Análisis de Sólidos Suspendidos

Autor/a

Afiliación

Ariadna Malena Seba

Fecha de publicación

8 de agosto de 2025

Resumen

Este sitio web aborda el análisis de sólidos suspendidos totales (SST) en cuerpos de agua mediante el uso de imágenes satelitales y técnicas de aprendizaje automático. Presenta una introducción teórica sobre la importancia de los SST como indicador ambiental y luego desarrolla una parte práctica con Python, donde se integran datos espectrales y mediciones reales para aplicar modelos de regresión. Se evalúa el desempeño de estos modelos y se exploran relaciones entre bandas espectrales y SST, proponiendo mejoras metodológicas para estudios futuros.

Palabras clave

GISTAQ, UTN, FRRe, Quarto, Sólidos Suspendidos

1 Sólidos suspendidos totales, por Vera Geneyer (https://github.com/VeraGeneyer)

Los sólidos suspendidos totales (TSM) son la cantidad de materia en suspensión en el agua, que incluye plancton, minerales, arena, y microorganismos. Se determinan como el residuo no filtrable de una muestra de agua. Niveles altos (TSM) pueden reducir la transparencia del agua, limitar la luz y y transportar sustancias tóxicas, afectando la vida acuática y la calidad del agua. Este parámetro, medido mediante sensores remotos, nos da información sobre el estado físico del cuerpo de agua y están relacionados con factores como la humedad, temperatura y entre otros, que es vital para detectar riesgos al ecosistema y cumplir con las normas ambientales.

1.1 Métodos tradicionales

Tabla 1: Características principales de algoritmos tradicionales para la estimación de sólidos suspendidos.

Ecuación	Bandas (nm)	Métricas	Aguas	Plataforma	Referencia
\(-229.34 \left( \frac{B03}{B08} \right)^{3}+1001.65 \left( \frac{B03}{B08} \right)^{2}-1422.7 \left( \frac{B03}{B08} \right)+665.17\)	B03, B08	\(R^{2}\)	Embalse¹	Landsat-8	[1]
\(-244.83+40.21 \cdot B01-3.67 \cdot NDWI\)	B01, NDWI (B03, B08)	\(R^{2}\), RMSE, d	Río²	GeoEye	[2]

De acuerdo a un estudio que analizó 48 cuerpos de agua, la estimación de TSM se hizo en su mayoría por modelos lineales, siendo la banda B8A la más frecuente [3].

1.2 Métodos de aprendizaje automático

El aprendizaje automático (ML) es una rama de la inteligencia artificial cuyo objetivo es desarrollar algoritmos capaces de resolver problemas mediante el análisis de datos y la creación de funciones que describen el comportamiento de fenómenos monitoreados [4]. Los modelos de aprendizaje automático más utilizados y mencionados por los investigadores para predecir la concentración de SST son:

Bosque Aleatorio (RF) y Refuerzo Adaptativo (AdB), modelos que se destacan por su robustez ante datos complejos y ruidosos. Estos algoritmos construyen múltiples árboles de decisión que analizan las relaciones entre características como el uso del suelo o el volumen de escorrentía y los niveles de SST [5].
Redes Neuronales Artificiales (ANN), copian las redes neuronales biológicas y aprenden patrones complejos en grandes volúmenes de datos, como los niveles de SST en distintas condiciones ambientales [5],
k-Nearest Neighbors (kNN), en sus variantes de ponderación uniforme y variable, que estima el SST en función de la cercanía en características de nuevos puntos de muestreo con datos históricos [5].

El aprendizaje automático es esencial para mejorar la precisión y rapidez en el análisis de la calidad del agua, proporcionando un monitoreo más eficiente y menos costoso en comparación con los métodos tradicionales, especialmente en áreas de difícil acceso o con datos limitados.

Tabla 2: Características principales de algoritmos de aprendizaje automático para la estimación de sólidos suspendidos.

Modelo de machine learning	Software	Agua	Datos	Métricas	Referencias
Bagging y Random Forest	Programa R	Bahía	Muestreo	Prueba de normalidad multivalente Mardia-tests y Royston	[4]
Regresión lineal, LASSO, regresión de vectores de soporte (SVR), K vecinos más cercanos (KNN), bosque aleatorio (RF) y redes neuronales artificiales (ANN).	-	Lago y embalse	Sentinel-2 y UAV	\(R^{2}\)	[6]
Regresión lineal, regresión de vectores de soporte (SVR), K vecinos más cercanos (KNN), bosque aleatorio (RF) y redes neuronales artificiales (ANN).	Programa Python	Lagos	Estación de monitoreo (Sensores para cada parámetro)	\(R^{2}\), NSE y RMSE	[5]

2 Desarrollo integral del algoritmo

En esta sección describo el proceso completo para crear un modelo que estima los sólidos suspendidos totales (SST) usando datos satelitales y de campo, desde la preparación de datos hasta la validación del algoritmo.

2.1 Procesamiento inicial de datos con Sen2Cor

Se detalla el procedimiento técnico que implementé para procesar información ambiental georreferenciada, con el objetivo de analizar el comportamiento del parámetro sólidos suspendidos (sol_sus) en una región específica del estudio (píxel 3x3).

Para esto, utilicé el lenguaje Python y la biblioteca pandas, que resulta particularmente eficiente para el manejo y transformación de datos en estructuras tabulares. Este paso fue fundamental para preparar los datos y asegurar su calidad antes del análisis.

2.1.1 Carga de datos

Primero importo la biblioteca pandas, una herramienta en Python que se utiliza para manejar datos en formato tabular (como hojas de cálculo o CSVs). Se le da el alias pd por convención, para simplificar el código.

Luego cargo dos archivos CSV con la función pd.read_csv(), la cual convierte dichos archivos en objetos del tipo DataFrame, que representan tablas en memoria, que son estructuras de datos similares a tablas (parecida a una hoja de Excel). Los conjuntos de datos cargados fueron:

gis_df: contiene información geográfica (latitud, longitud, pixel, etc.).
lab_df: contiene datos de laboratorio, incluyendo el parámetro de interés sol_sus.

Verifico la carga correcta mostrando las primeras filas con la función .head(). Es útil para ver rápidamente cómo es la estructura del archivo: qué columnas hay, qué tipo de datos, si se cargó bien.

import pandas as pd  # pandas es la biblioteca para manejar datos tabulares

# Cargar los archivos de datos
gis_df = pd.read_csv('datos/base_de_datos_gis.csv')
lab_df = pd.read_csv('datos/base_de_datos_lab.csv')

# Ver las primeras filas para asegurarse de que se cargaron bien
gis_df.head(), lab_df.head()

print("Primeras filas de gis_df:")
display(gis_df.head())

print("\nPrimeras filas de lab_df:")
display(lab_df.head())

Primeras filas de gis_df:

	fecha	punto	pixel	banda	reflect	longitud	latitud
0	2023-05-11	1	1x1	B01	0.161744	-58.868047	-27.464687
1	2023-05-11	1	1x1	B02	0.170400	-58.868047	-27.464687
2	2023-05-11	1	1x1	B03	0.208800	-58.868047	-27.464687
3	2023-05-11	1	1x1	B04	0.251000	-58.868047	-27.464687
4	2023-05-11	1	1x1	B05	0.257656	-58.868047	-27.464687


Primeras filas de lab_df:

	fecha	longitud	latitud	param	valor
0	2023-05-11	-58.868047	-27.464687	ph	6.8
1	2023-05-11	-58.868047	-27.464687	cond	141.1
2	2023-05-11	-58.868047	-27.464687	sol_sus	198.0
3	2023-05-11	-58.868047	-27.464687	turb	185.5
4	2023-05-11	-58.866729	-27.466303	ph	6.8

2.1.2 Filtrado del parámetro sólidos suspendidos (‘sol_sus’)

En el conjunto de datos del laboratorio lab_df, hay múltiples parámetros medidos (como pH, turbidez, etc.). En este caso, me interesa trabajar solamente con los datos de sólidos suspendidos, identificado como "sol_sus" en la columna param. Este filtrado selectivo lo realicé para limitar el análisis al fenómeno físico-químico de interés.

Filtré el DataFrame para quedarme solo con esas filas, y renombré la columna valor como sol_sus para que sea más claro en los siguientes pasos.

# Filtrar solo las filas donde el parámetro es 'sol_sus'
sol_sus_df = lab_df[lab_df["param"] == "sol_sus"]

# Renombrar la columna 'valor' a 'sol_sus' para que tenga sentido en el merge
sol_sus_df = sol_sus_df.rename(columns={"valor": "sol_sus"})

# Mostrar para confirmar
sol_sus_df.head()
print("Primeras filas de sol_sus_df:")
display(sol_sus_df.head())

Primeras filas de sol_sus_df:

	fecha	longitud	latitud	param	sol_sus
2	2023-05-11	-58.868047	-27.464687	sol_sus	198.0
6	2023-05-11	-58.866729	-27.466303	sol_sus	150.0
10	2023-05-11	-58.864889	-27.468056	sol_sus	101.0
14	2023-05-11	-58.863367	-27.469240	sol_sus	95.0
18	2023-05-11	-58.862454	-27.470561	sol_sus	69.0

2.1.3 Reorganización de datos: pivotar bandas en columnas individuales

En este paso, convierto los valores únicos de la columna banda (como B01, B02, etc.) en nombres de columnas. Cada nueva columna contendrá los valores del parámetro reflect correspondientes a esa banda en particular. Esta operación se realiza antes de unir con los valores de sol_sus, ya que el valor de reflectancia depende de la banda, mientras que sol_sus es un dato independiente que se asignará luego por punto, fecha y ubicación.

# Pivotear la tabla para que cada banda sea una columna
gis_pivot = gis_df.pivot_table(
    index=['fecha', 'punto', 'pixel', 'latitud', 'longitud'],
    columns='banda',
    values='reflect'
).reset_index()

# Eliminar el nombre del índice de columnas generado por el pivot
gis_pivot.columns.name = None

print("Primeras filas de gis_pivot:")
display(gis_pivot.head())

Primeras filas de gis_pivot:

	fecha	punto	pixel	latitud	longitud	B01	B02	B03	B04	B05	B06	B07	B08	B11	B12	B8A
0	2023-05-11	1	1x1	-27.464687	-58.868047	0.161744	0.170400	0.208800	0.251000	0.257656	0.202175	0.207500	0.194000	0.118625	0.114738	0.173069
1	2023-05-11	1	3x3	-27.464687	-58.868047	0.161290	0.171256	0.206644	0.250911	0.258058	0.202170	0.207490	0.193722	0.118892	0.115462	0.173137
2	2023-05-11	2	1x1	-27.466303	-58.866729	0.149925	0.170600	0.202100	0.243200	0.247169	0.186631	0.187512	0.175800	0.114100	0.111031	0.159263
3	2023-05-11	2	3x3	-27.466303	-58.866729	0.149906	0.170378	0.202278	0.242689	0.247115	0.186739	0.188027	0.175544	0.113992	0.111210	0.159373
4	2023-05-11	3	1x1	-27.468056	-58.864889	0.149725	0.166200	0.197100	0.223800	0.218175	0.159169	0.163344	0.150700	0.114263	0.113006	0.139662

2.1.4 Unión de datos geoespaciales con datos de laboratorio

Una vez que las bandas han sido transformadas en columnas, combino esta tabla con los valores de sólidos suspendidos (sol_sus) provenientes del laboratorio. La combinación se hace usando las columnas fecha, latitud y longitud, que permiten identificar los datos correspondientes a un mismo punto geográfico y temporal.

# Realizar el merge por ubicación y fecha
df_merged = pd.merge(
    gis_pivot,
    sol_sus_df[['fecha', 'latitud', 'longitud', 'sol_sus']],
    on=['fecha', 'latitud', 'longitud'],
    how='inner'
)

print("Primeras filas del DataFrame combinado:")
display(df_merged.head())

Primeras filas del DataFrame combinado:

	fecha	punto	pixel	latitud	longitud	B01	B02	B03	B04	B05	B06	B07	B08	B11	B12	B8A	sol_sus
0	2023-05-11	1	1x1	-27.464687	-58.868047	0.161744	0.170400	0.208800	0.251000	0.257656	0.202175	0.207500	0.194000	0.118625	0.114738	0.173069	198.0
1	2023-05-11	1	3x3	-27.464687	-58.868047	0.161290	0.171256	0.206644	0.250911	0.258058	0.202170	0.207490	0.193722	0.118892	0.115462	0.173137	198.0
2	2023-05-11	2	1x1	-27.466303	-58.866729	0.149925	0.170600	0.202100	0.243200	0.247169	0.186631	0.187512	0.175800	0.114100	0.111031	0.159263	150.0
3	2023-05-11	2	3x3	-27.466303	-58.866729	0.149906	0.170378	0.202278	0.242689	0.247115	0.186739	0.188027	0.175544	0.113992	0.111210	0.159373	150.0
4	2023-05-11	3	1x1	-27.468056	-58.864889	0.149725	0.166200	0.197100	0.223800	0.218175	0.159169	0.163344	0.150700	0.114263	0.113006	0.139662	101.0

2.1.5 Filtrado espacial por píxel de interés (3x3)

Luego de combinar los datos, aplico un filtrado adicional al DataFrame sobre la columna pixel para conservar únicamente las filas correspondientes al área geográfica designada como "3x3". Este paso reduce el dominio de análisis y permite concentrarse en una región de estudio concreta.

# Filtrar solo los datos del pixel 3x3
df_pixel_3x3 = df_merged[df_merged["pixel"] == "3x3"]

print("Primeras filas del pixel 3x3:")
display(df_pixel_3x3.head())

Primeras filas del pixel 3x3:

	fecha	punto	pixel	latitud	longitud	B01	B02	B03	B04	B05	B06	B07	B08	B11	B12	B8A	sol_sus
1	2023-05-11	1	3x3	-27.464687	-58.868047	0.161290	0.171256	0.206644	0.250911	0.258058	0.202170	0.207490	0.193722	0.118892	0.115462	0.173137	198.0
3	2023-05-11	2	3x3	-27.466303	-58.866729	0.149906	0.170378	0.202278	0.242689	0.247115	0.186739	0.188027	0.175544	0.113992	0.111210	0.159373	150.0
5	2023-05-11	3	3x3	-27.468056	-58.864889	0.149718	0.167256	0.197322	0.223511	0.218608	0.159292	0.163228	0.151544	0.114241	0.112931	0.139783	101.0
7	2023-05-11	4	3x3	-27.469240	-58.863367	0.148441	0.163289	0.188733	0.196167	0.183854	0.139388	0.140017	0.135400	0.115602	0.113097	0.128222	95.0
9	2023-05-11	5	3x3	-27.470561	-58.862454	0.145838	0.157733	0.178767	0.173911	0.162600	0.132140	0.133229	0.126656	0.117069	0.113992	0.124548	69.0

2.1.6 Exportación de datos procesados

Finalmente, guardo el resultado como un nuevo archivo .csv dentro de la carpeta datos.

Por último, exporto el resultado a un nuevo archivo en formato .csv, mediante la función to_csv() de pandas, con el parámetro index=False para evitar que la columna de índice se incluya en el archivo de salida que pandas crea por defecto. Esto me permite utilizarlo después para visualización o análisis posterior.

# Guardar el archivo CSV dentro de la carpeta "datos"
df_pixel_3x3.to_csv('datos/datos_sol_sus_pixel_3x3.csv', index=False)

2.2 Aplicación del procesamiento a datos ACOLITE

A continuación, repetí el mismo procedimiento de procesamiento aplicado anteriormente, esta vez utilizando el archivo base_de_datos_gis_acolite.csv, generado con el procesador atmosférico ACOLITE.

ACOLITE es una herramienta diseñada específicamente para corregir efectos atmosféricos en ambientes acuáticos, a partir de imágenes satelitales. Gracias a que la estructura del archivo es similar al utilizado anteriormente (Sen2Cor), fue posible aplicar la misma lógica de filtrado, transformación y combinación para preparar los datos de reflectancia.

Código

import pandas as pd

# 1. Cargar archivos CSV 
gis_acol_df = pd.read_csv('datos/base_de_datos_gis_acolite.csv')
lab_df = pd.read_csv('datos/base_de_datos_lab.csv')

# 2. Filtrar solo el parámetro 'sol_sus' 
sol_sus_df = lab_df[lab_df["param"] == "sol_sus"]
sol_sus_df = sol_sus_df.rename(columns={"valor": "sol_sus"})

#  3. Pivotear bandas en columnas 
gis_pivot = gis_acol_df.pivot_table(
    index=['fecha', 'punto', 'latitud', 'longitud'],
    columns='banda',
    values='reflect'
).reset_index()

gis_pivot.columns.name = None  # Quitar la etiqueta 'banda'

# 4. Combinar reflectancias con sol_sus 
df_merged = pd.merge(
    gis_pivot,
    sol_sus_df[['fecha', 'latitud', 'longitud', 'sol_sus']],
    on=['fecha', 'latitud', 'longitud'],
    how='inner'
)

# 5. Guardar resultado final
df_merged.to_csv('datos/datos_sol_sus_acolite.csv', index=False)

# 6. Mostrar muestra del resultado
print("Primeras filas del DataFrame resultante:")
display(df_merged.head())

Primeras filas del DataFrame resultante:

	fecha	punto	latitud	longitud	B01	B02	B03	B04	B05	B06	B07	B08	B11	B12	B8A	sol_sus
0	2023-05-11	1	-27.464687	-58.868047	0.041367	0.046669	0.085630	0.128799	0.128704	0.091748	0.098365	0.083253	0.002521	0.000000	0.066612	198.0
1	2023-05-11	2	-27.466303	-58.866729	0.032235	0.047916	0.084191	0.123898	0.122481	0.079922	0.082403	0.069021	0.001669	0.000000	0.055924	150.0
2	2023-05-11	3	-27.468056	-58.864889	0.030956	0.044084	0.078444	0.104498	0.095737	0.052162	0.057475	0.044354	0.000209	0.000174	0.034767	101.0
3	2023-05-11	4	-27.469240	-58.863367	0.029646	0.040407	0.070384	0.078357	0.063937	0.032452	0.033994	0.028441	0.001279	0.000000	0.022597	95.0
4	2023-05-11	5	-27.470561	-58.862454	0.026737	0.034968	0.060701	0.056602	0.044067	0.024809	0.026684	0.019209	0.001879	0.000000	0.017930	69.0

2.3 Modelo base · Regresión lineal simple

En este análisis aplico un modelo de regresión lineal simple para estudiar la relación entre la reflectancia (banda B01) y la concentración de sólidos suspendidos (sol_sus), utilizando datos experimentales. La regresión lineal es una técnica fundamental del aprendizaje automático supervisado que permite predecir un valor continuo a partir de una o más variables independientes, y suele emplearse como modelo base frente a enfoques más complejos. A lo largo de esta sección, detallo paso a paso las acciones realizadas y explico los conceptos clave para comprender y replicar este análisis.

2.3.1 Importar librerías

En este paso, cargo las bibliotecas necesarias para procesar datos, ajustar modelos de regresión, evaluar su desempeño y visualizar los resultados.

import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt

2.3.2 Cargar datos desde un CSV

Importo el archivo .csv con los datos experimentales. Se visualizan las primeras filas para verificar que los datos se han cargado correctamente.

# Cargar el CSV
datos = pd.read_csv('datos/datos_sol_sus_acolite.csv')

# Mostrar las primeras filas para verificar
datos.head()
print("Primeras filas de datos:")
display(datos.head())

Primeras filas de datos:

	fecha	punto	latitud	longitud	B01	B02	B03	B04	B05	B06	B07	B08	B11	B12	B8A	sol_sus
0	2023-05-11	1	-27.464687	-58.868047	0.041367	0.046669	0.085630	0.128799	0.128704	0.091748	0.098365	0.083253	0.002521	0.000000	0.066612	198.0
1	2023-05-11	2	-27.466303	-58.866729	0.032235	0.047916	0.084191	0.123898	0.122481	0.079922	0.082403	0.069021	0.001669	0.000000	0.055924	150.0
2	2023-05-11	3	-27.468056	-58.864889	0.030956	0.044084	0.078444	0.104498	0.095737	0.052162	0.057475	0.044354	0.000209	0.000174	0.034767	101.0
3	2023-05-11	4	-27.469240	-58.863367	0.029646	0.040407	0.070384	0.078357	0.063937	0.032452	0.033994	0.028441	0.001279	0.000000	0.022597	95.0
4	2023-05-11	5	-27.470561	-58.862454	0.026737	0.034968	0.060701	0.056602	0.044067	0.024809	0.026684	0.019209	0.001879	0.000000	0.017930	69.0

2.3.3 Definir variables y particionar el conjunto

Selecciono las variables relevantes: B01 como variable independiente y sol_sus como variable dependiente. Luego divido el conjunto en dos subconjuntos: uno para entrenar el modelo y otro para probarlo, lo cual sirve para evaluar su capacidad de generalización.

# Selección de variables
X = datos[["B01"]]           # variable independiente
y = datos["sol_sus"]         # variable dependiente

# División en entrenamiento y validación
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=False)

2.3.4 Entrenar modelo de regresión lineal

En este paso se entrena un modelo de regresión lineal usando los datos de entrenamiento. El modelo aprende la relación matemática entre la reflectancia y los sólidos suspendidos.

regressor = LinearRegression().fit(X_train, y_train)

2.3.5 Evaluar desempeño del modelo

Una vez entrenado el modelo, evaluo su desempeño usando métricas estadísticas. Estas nos permiten cuantificar qué tan bien el modelo predice los valores de sólidos suspendidos a partir de la reflectancia en los datos de prueba.

y_pred = regressor.predict(X_test)
p_rmse = mean_squared_error(y_test, y_pred)
p_r2 = r2_score(y_test, y_pred)

Métricas de desempeño

El error cuadrático medio es: 1846.79
El coeficiente de determinación (R²) es: -0.252

2.3.6 Visualizar resultados

Finalmente, se visualiza gráficamente la relación entre reflectancia y sólidos suspendidos, tanto en el conjunto de entrenamiento como en el de prueba. Esto ayuda a interpretar de forma visual cómo se ajusta el modelo a los datos reales.

fig, ax = plt.subplots(ncols=2, figsize=(10, 5), sharex=True, sharey=True)

# Gráfico entrenamiento
ax[0].plot(X_train, regressor.predict(X_train), linewidth=3, color="#17A77E", label="Modelo")
ax[0].scatter(X_train, y_train, label="Entrenamiento", color="#9D50A6", alpha=0.6)
ax[0].set(xlabel="Reflectancia", ylabel="Sol_Sus", title="Conjunto de entrenamiento")
ax[0].legend()

# Gráfico validación
ax[1].plot(X_test, y_pred, linewidth=3, color="#17A77E", label="Modelo")
ax[1].scatter(X_test, y_test, label="Validación", color="#9D50A6", alpha=0.6)
ax[1].set(xlabel="Reflectancia", ylabel="Sol_Sus", title="Conjunto de validación")
ax[1].legend()

fig.suptitle("Regresión lineal")

plt.show()

2.4 Regresión lineal individual por banda

Con el fin de profundizar el análisis, evalué la relación entre los sólidos suspendidos (sol_sus) y cada banda espectral por separado. Para esto, entrené un modelo de regresión lineal simple usando los mismos datos experimentales para cada banda.

Con este enfoque comparo la capacidad predictiva individual de cada banda mediante las métricas R², R² ajustado y el error cuadrático medio (RMSE). Este análisis me permite identificar qué bandas tienen mejor relación lineal con los sólidos suspendidos y sentar la base para probar modelos más complejos.

Código

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
import math
from IPython.display import display

# Cargar datos
datos = pd.read_csv('datos/datos_sol_sus_acolite.csv')

# Detectar columnas de bandas
bandas = [col for col in datos.columns if col.startswith("B")]

# Lista para guardar resultados
resultados = []

# Parámetros para organización de gráficos
n_bandas = len(bandas)
ncols = 3  # Número de columnas de la grilla
nrows = math.ceil(n_bandas / ncols)  # Calculamos cuántas filas se necesitan

# Crear figura
fig, axs = plt.subplots(nrows=nrows, ncols=ncols, figsize=(12, 4 * nrows))
axs = axs.flatten()  # Asegura que podamos indexarlos como una lista

for i, banda in enumerate(bandas):
    # Variables
    X = datos[[banda]]
    y = datos["sol_sus"]

    # División entrenamiento/prueba
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=False)

    # Ajuste del modelo
    modelo = LinearRegression().fit(X_train, y_train)
    y_train_pred = modelo.predict(X_train)

    # Métricas
    r2 = modelo.score(X_train, y_train)
    n = len(y_train)
    p = X_train.shape[1]
    r2_adj = 1 - (1 - r2) * (n - 1) / (n - p - 1)
    rmse = np.sqrt(mean_squared_error(y_train, y_train_pred))

    resultados.append({
        "Banda": banda,
        "R²": round(r2, 4),
        "R²_ajustado": round(r2_adj, 4),
        "RMSE": round(rmse, 4)
    })

    # Gráfico de entrenamiento
    ax = axs[i]
    ax.scatter(X_train, y_train, alpha=0.6, color="#9D50A6", label="Entrenamiento")
    ax.plot(X_train, y_train_pred, color="#17A77E", linewidth=1.8, label="Modelo")
    ax.set_title(f'{banda}\nR²={r2:.2f}', fontsize=10)
    ax.set_xlabel('Reflectancia', fontsize=8)
    ax.set_ylabel('Sol_Sus', fontsize=8)
    ax.tick_params(axis='both', which='major', labelsize=8)
    ax.legend(fontsize=7)
    ax.grid(True)

# Eliminar ejes sobrantes
for j in range(i + 1, len(axs)):
    fig.delaxes(axs[j])

plt.suptitle("Regresiones lineales por banda (entrenamiento)", fontsize=14)
plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

# Tabla resumen
df_resultados = (
    pd.DataFrame(resultados)
      .round(3)
      .sort_values("R²", ascending=False)   # la mejor banda queda arriba
      .reset_index(drop=True)
)

def style_metrics(df: pd.DataFrame, caption: str):
    sty = (
        df.style
          .format(precision=3)
    )

    # Ocultar índice según versión de pandas
    if hasattr(sty, "hide"):
        sty = sty.hide(axis="index")
    else:
        sty = sty.hide_index()

    # Agregar título y estilo del título
    sty = sty.set_caption(caption)
    sty = sty.set_table_styles([{
        "selector": "caption",
        "props": [
            ("caption-side", "top"),
            ("font-weight",  "bold"),
            ("margin-bottom", "0.5em")
        ]
    }], overwrite=False)

    return sty

# Mostrar tabla  
display(style_metrics(df_resultados,
                      "Resumen de métricas por banda (entrenamiento)"))

Tabla 3: Resumen de métricas por banda (entrenamiento)

Banda	R²	R²_ajustado	RMSE
B05	0.183	0.167	33.457
B06	0.180	0.164	33.523
B07	0.173	0.157	33.651
B08	0.146	0.130	34.191
B8A	0.135	0.118	34.429
B04	0.094	0.076	35.230
B12	0.038	0.019	36.303
B11	0.033	0.015	36.389
B01	0.008	-0.011	36.865
B03	0.005	-0.014	36.917
B02	0.001	-0.018	36.985

Conclusiones del análisis por banda

Los resultados muestran que ninguna banda espectral logra predecir con alta precisión los sólidos suspendidos mediante un modelo lineal simple.

Las bandas B05, B06 y B07 presentan los mejores valores de R² (alrededor de 0.18), lo que indica que explican menos del 20 % de la variabilidad en los datos. El resto de las bandas tienen un desempeño aún menor, con errores (RMSE) relativamente altos en comparación con los valores observados.

Esto sugiere que la relación entre reflectancia y sólidos suspendidos no es lineal en escala natural y que será necesario explorar modelos más complejos o transformaciones de datos para mejorar la predicción.

2.5 Análisis de regresión lineal con transformación logarítmica

En esta etapa aplico una transformación logarítmica natural a las variables de reflectancia y sólidos suspendidos antes de ajustar los modelos de regresión lineal. Esta transformación permite estabilizar la varianza, linealizar relaciones no lineales y reducir el impacto de valores extremos.

El procedimiento es similar al análisis anterior, pero antes de entrenar el modelo se aplica log(x) a las columnas correspondientes. Para evitar problemas con valores cero, estos se reemplazan por NaN y se excluyen del análisis. Posteriormente, entreno modelos de regresión lineal simple por banda usando las variables transformadas.

El desempeño de cada modelo se evalúa mediante métricas en escala logarítmica (R²_log, RMSE_log) y en escala original, facilitando la comparación y la interpretación de resultados.

Código

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import math
import seaborn as sns
from IPython.display import display

# 1. ── Cargar datos ────────────────────────────────────────────────
datos = pd.read_csv("datos/datos_sol_sus_acolite.csv")
bandas = [c for c in datos.columns if c.startswith("B")]

# 2. ── Filtrar positivos (requisito para log) ──────────────────────
mask_pos   = (datos["sol_sus"] > 0) & (datos[bandas] > 0).all(axis=1)
datos_log  = datos.loc[mask_pos].copy()

# 3. ── Transformar a log-natural ──────────────────────────────────
datos_log["log_sol_sus"] = np.log(datos_log["sol_sus"])
for b in bandas:
    datos_log[f"log_{b}"] = np.log(datos_log[b])

log_bandas = [f"log_{b}" for b in bandas]

# 4. ── Train / test split (descartamos test) ───────────────────────
X_total = datos_log[log_bandas]
y_total = datos_log["log_sol_sus"]
X_train, _, y_train, _ = train_test_split(
    X_total, y_total, test_size=0.2, shuffle=False)

# 5. ── Funciones de métricas ───────────────────────────────────────
def metrics_log(model, X, y):
    y_pred = model.predict(X)
    rmse   = np.sqrt(mean_squared_error(y, y_pred))
    r2     = r2_score(y, y_pred)
    n, p   = X.shape
    r2_adj = 1 - ((1 - r2) * (n - 1)) / (n - p - 1)
    return rmse, r2, r2_adj

def metrics_original(model, X, y):
    y_pred_log  = model.predict(X).ravel()
    y_pred_orig = np.exp(y_pred_log)
    y_true_orig = np.exp(y.values)
    rmse   = np.sqrt(mean_squared_error(y_true_orig, y_pred_orig))
    r2     = r2_score(y_true_orig, y_pred_orig)
    n, p   = X.shape
    r2_adj = 1 - ((1 - r2) * (n - 1)) / (n - p - 1)
    return rmse, r2, r2_adj

# 6. ── Entrenar un modelo por banda ────────────────────────────────
resultados, resultados_nolog = [], []
lineas_log, lineas_nolog     = {}, {}

for b, lb in zip(bandas, log_bandas):
    Xtr = X_train[[lb]]
    reg = LinearRegression().fit(Xtr, y_train)

    # Métricas en log
    rmse_l, r2_l, r2a_l = metrics_log(reg, Xtr, y_train)
    resultados.append({"Banda": b, "R²_log": r2_l, 
                    "R²aj_log": r2a_l, "RMSE_log": rmse_l})

    # Métricas en original
    rmse_o, r2_o, r2a_o = metrics_original(reg, Xtr, y_train)
    resultados_nolog.append({"Banda": b,
                             "R²": r2_o, "R²aj": r2a_o, "RMSE": rmse_o})

    # Rectas ordenadas
    order        = np.argsort(Xtr.values.ravel())
    x_log_ord    = Xtr.values.ravel()[order]
    y_pred_log   = reg.predict(Xtr).ravel()[order]
    lineas_log[b] = (x_log_ord, y_pred_log)

    x_orig_ord   = np.exp(x_log_ord)
    y_pred_orig  = np.exp(y_pred_log)
    lineas_nolog[b] = (x_orig_ord, y_pred_orig)

# 7. ── DataFrames de métricas ──────────────────────────────────────
res_df_log   = (pd.DataFrame(resultados)
                  .round(3)
                  .sort_values("RMSE_log")
                  .reset_index(drop=True))

res_df_nolog = (pd.DataFrame(resultados_nolog)
                  .round(3)
                  .sort_values("RMSE")
                  .reset_index(drop=True))

# 8. ── Gráficos en escala original ────────────────────────────────
n_cols, n_rows = 3, math.ceil(len(bandas) / 3)
fig, axes = plt.subplots(n_rows, n_cols, figsize=(15, 4 * n_rows))
axes = axes.flatten()

for idx, b in enumerate(bandas):
    ax = axes[idx]
    x_orig = np.exp(X_train[f"log_{b}"])
    y_orig = np.exp(y_train)
    ax.scatter(x_orig, y_orig, alpha=0.6, color="#9D50A6", label="Entrenamiento")
    ax.plot(*lineas_nolog[b], lw=3, color="#17A77E", label="Modelo")
    ax.set(xlabel=b, ylabel="sol_sus", title=f"Banda {b} (escala original)")
    ax.legend()
    ax.grid(True)

for j in range(len(bandas), len(axes)):
    fig.delaxes(axes[j])

plt.tight_layout()
plt.show()

# 9. ── Función para estilo de tablas ───────────────────────────────
def style_metrics_sin_gradiente(df: pd.DataFrame, caption: str):
    sty = (df.style
             .format(precision=3)
             # .background_gradient(...) ← quitamos el gradiente
             .hide(axis="index")
             .set_caption(caption)
             .set_table_styles([{
                 "selector": "caption",
                 "props": [("caption-side", "top"),
                           ("font-weight",  "bold"),
                           ("margin-bottom","0.5em")]
             }], overwrite=False)
          )
    return sty


# 10 ── Mostrar tablas con el nuevo estilo pastel ───────────────────────────────
display(style_metrics_sin_gradiente(res_df_log,
        "Métricas en escala log‑log (ordenadas por RMSE_log)"))

display(style_metrics_sin_gradiente(res_df_nolog,
        "Métricas en escala original (ordenadas por RMSE)"))

Tabla 4: Métricas en escala log‑log (ordenadas por RMSE_log)

Banda	R²_log	R²aj_log	RMSE_log
B06	0.304	0.281	0.318
B07	0.303	0.280	0.319
B05	0.272	0.248	0.326
B08	0.267	0.243	0.327
B8A	0.262	0.237	0.328
B04	0.166	0.138	0.349
B12	0.111	0.082	0.360
B11	0.065	0.034	0.369
B03	0.017	-0.016	0.378
B01	0.015	-0.018	0.379
B02	0.007	-0.026	0.380

Tabla 5: Métricas en escala original (ordenadas por RMSE)

Banda	R²	R²aj	RMSE
B06	0.198	0.171	31.661
B07	0.186	0.159	31.897
B05	0.178	0.150	32.049
B8A	0.156	0.128	32.472
B08	0.151	0.123	32.567
B04	0.077	0.046	33.962
B12	0.076	0.045	33.986
B11	0.043	0.011	34.576
B01	-0.028	-0.063	35.843
B03	-0.034	-0.068	35.940
B02	-0.035	-0.070	35.969

Conclusión del ajuste log–log por banda

Los resultados muestran que al aplicar la transformación logarítmica, las bandas B06, B07 y B05 obtienen los mejores ajustes, con valores de R²_log cercanos a 0.30 y errores (RMSE_log) más bajos comparados con las demás bandas. Esto indica una mejora en la capacidad predictiva respecto al análisis en escala natural.

Sin embargo, al evaluar las métricas en escala original, el mejor desempeño alcanza un R² máximo cercano a 0.20, lo que aún refleja una capacidad limitada para explicar la variabilidad de los sólidos suspendidos utilizando bandas individuales.

Estos resultados sugieren que la transformación logarítmica mejora el modelo, pero que para obtener predicciones más precisas será necesario explorar combinaciones de bandas, índices espectrales o modelos más complejos.

2.6 Selección de bandas mediante regresión logarítmica y AIC

Para determinar cuáles bandas espectrales aportan más información para predecir sólidos suspendidos, ajusté modelos de regresión lineal simple en escala logarítmica para cada banda individual. Además de las métricas habituales como el coeficiente de determinación (R²) y el error cuadrático medio (RMSE), utilicé el Criterio de Información de Akaike (AIC).

El AIC es una medida que evalúa la calidad del modelo penalizando la complejidad: un valor menor indica un mejor equilibrio entre ajuste y simplicidad. Esto ayuda a evitar sobreajuste al no favorecer modelos que mejoran el ajuste sólo por aumentar el número de parámetros sin aportar realmente información útil. Así, el AIC es una herramienta clave para seleccionar las bandas que efectivamente mejoran la predicción sin agregar complejidad innecesaria.

Código

import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error, r2_score
import seaborn as sns
from IPython.display import display

# 1. Cargar datos y detectar bandas
datos = pd.read_csv("datos/datos_sol_sus_acolite.csv")
bandas = [c for c in datos.columns if c.startswith("B")]

# 2. Filtrar casos con valores positivos (requisito de log)
mask_pos = (datos["sol_sus"] > 0) & (datos[bandas] > 0).all(axis=1)
datos_log = datos.loc[mask_pos].copy()

# 3. Transformar a log-natural
datos_log["log_sol_sus"] = np.log(datos_log["sol_sus"])
for b in bandas:
    datos_log[f"log_{b}"] = np.log(datos_log[b])

log_bandas = [f"log_{b}" for b in bandas]

# 4. Partición entrenamiento / test 
X_total = datos_log[log_bandas]
y_total = datos_log["log_sol_sus"]

X_train, _, y_train, _ = train_test_split(
    X_total, y_total, test_size=0.20, shuffle=False
)

# 5. Funciones auxiliares (AIC + métricas)
def aic_gauss(model, X, y):
    """
    AIC para regresión lineal con ruido ~ N(0, σ²)
    AIC = n · ln(RSS/n) + 2k
    """
    n = len(y)
    k = X.shape[1] + 1  # predictores + intercepto
    rss = np.sum((y - model.predict(X)) ** 2)
    return n * np.log(rss / n) + 2 * k

def metrics_log(model, X, y):
    y_pred = model.predict(X)
    rmse = np.sqrt(mean_squared_error(y, y_pred))
    r2   = r2_score(y, y_pred)
    n, p = X.shape
    r2_adj = 1 - ((1 - r2) * (n - 1)) / (n - p - 1)
    return rmse, r2, r2_adj

# 6. Ajuste univariable banda-a-banda + AIC
resultados = []

for b in bandas:
    lb = f"log_{b}"
    Xtr = X_train[[lb]]
    reg = LinearRegression().fit(Xtr, y_train)

    rmse_log, r2_log, r2a_log = metrics_log(reg, Xtr, y_train)
    aic_val = aic_gauss(reg, Xtr, y_train)

    resultados.append({
        "Banda"    : b,
        "R²"   : r2_log,
        "R²aj" : r2a_log,
        "RMSE" : rmse_log,
        "AIC"      : aic_val
    })

# 7. Tabla ordenada por AIC
df_log  = (pd.DataFrame(resultados)
             .round(3)
             .sort_values("AIC")
             .reset_index(drop=True))

def style_metrics_sin_gradiente(df: pd.DataFrame, caption: str):
    sty = (df.style
             .format(precision=3)
             .hide(axis="index")  # oculta índice
             .set_caption(caption)
             .set_table_styles([{
                 "selector": "caption",
                 "props": [("caption-side", "top"),
                           ("font-weight",  "bold"),
                           ("margin-bottom", "0.5em")]
             }], overwrite=False)
          )
    return sty

# 9 ── Mostrar tabla con el nuevo estilo pastel AIC ───────────────────────────
display(style_metrics_sin_gradiente(
    df_log, "AIC y métricas en escala log‑log (ordenadas por AIC)"
))

Tabla 6: AIC y métricas en escala log‑log (ordenadas por AIC)

Banda	R²	R²aj	RMSE	AIC
B06	0.304	0.281	0.318	-69.242
B07	0.303	0.280	0.319	-69.191
B05	0.272	0.248	0.326	-67.798
B08	0.267	0.243	0.327	-67.584
B8A	0.262	0.237	0.328	-67.355
B04	0.166	0.138	0.349	-63.448
B12	0.111	0.082	0.360	-61.415
B11	0.065	0.034	0.369	-59.781
B03	0.017	-0.016	0.378	-58.184
B01	0.015	-0.018	0.379	-58.112
B02	0.007	-0.026	0.380	-57.878

Conclusión sobre la selección de bandas con AIC

Los valores obtenidos muestran que las bandas B06 y B07 son las más relevantes para predecir los sólidos suspendidos, ya que presentan los valores más bajos de AIC (-69.24 y -69.19 respectivamente), indicando el mejor equilibrio entre calidad de ajuste y simplicidad del modelo.
Estas bandas también exhiben los mejores coeficientes de determinación (R²_log alrededor de 0.30) y los menores errores (RMSE_log cerca de 0.32).

El resto de las bandas tienen valores de AIC menos favorables, lo que sugiere que, individualmente, aportan menos información útil y no justifican la complejidad del modelo. Esto refuerza la idea de que para mejorar las predicciones será necesario combinar varias bandas o explorar modelos más complejos.

2.6.1 Selección de banda espectral con validación cruzada

En esta sección selecciono la banda espectral más adecuada para estimar los sólidos suspendidos (sol_sus) usando un modelo de regresión lineal simple en escala logarítmica. A diferencia de análisis anteriores, ahora incorporo una validación cruzada que respeta la estructura temporal de los datos, para asegurar que el modelo sea confiable al aplicarse a nuevas campañas de muestreo.

¿Por qué validación Leave-One-Group-Out (LOGO)?

Los datos provienen de campañas realizadas en diferentes fechas (fecha). Esto genera grupos independientes, ya que las condiciones del río pueden variar entre campañas. Usar una validación cruzada tradicional como K-Fold o Shuffle Split mezclaría muestras de una misma campaña entre entrenamiento y prueba, lo que daría una estimación demasiado optimista del desempeño.

En cambio, LeaveOneGroupOut (LOGO) deja una fecha completa afuera en cada iteración, lo que permite evaluar si el modelo puede generalizar a campañas no vistas. Esta estrategia evita filtraciones de información (data leakage) entre los datos de entrenamiento y validación, y resulta mucho más realista en contextos ambientales, donde el objetivo es predecir nuevas campañas a partir de mediciones anteriores.

Descripción del procedimiento:

Preprocesamiento: aplico log() a las bandas y a la variable objetivo (sol_sus) y elimino ceros, para linearizar y estabilizar la relación.
Modelado: ajusto un modelo de regresión lineal simple por cada banda, normalizando los datos con StandardScaler (pipeline). Para cada modelo calculo las métricas de validación cruzada: RMSE, R², R² ajustado y el AIC, que penaliza la complejidad del modelo.
Selección: ordeno las bandas por AIC y RMSE, y elijo la mejor.
Evaluación final:
- Estimo el desempeño real sobre la campaña más reciente (test).
- Realizo un análisis bootstrap (500 réplicas) para obtener intervalos de confianza de los parámetros y métricas.
- Muestro gráficos de predicho vs observado en escala log y original.

Este enfoque no solo me permite elegir la banda más informativa, sino que también evalúa la estabilidad del modelo y su capacidad de generalizar en el tiempo, algo clave para el monitoreo operativo de calidad de agua.

Código

# -------------------------------------------------------------
# 0 · LIBRERÍAS
import numpy as np
import pandas as pd
from pathlib import Path
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import StandardScaler
from sklearn.pipeline import Pipeline
from sklearn.model_selection import LeaveOneGroupOut
from sklearn.metrics import mean_squared_error, r2_score
from IPython.display import display, Markdown

# ----------  ESTÉTICA TABLAS (sin gradiente) ------------------
def style_table(df, caption):
    return (df.style
              .format(precision=3)
              .hide(axis="index")
              .set_caption(caption)
              .set_table_styles([
                 {"selector": "caption",
                  "props": [("caption-side", "top"),
                            ("font-weight",  "bold"),
                            ("margin-bottom","0.5em"),
                            ("text-align","left")]},
                 {"selector": "th",
                  "props": [("text-align","center")]},
                 {"selector": "td",
                  "props": [("text-align","right"),
                            ("padding","0.25em 0.5em")]}
              ], overwrite=False))

# -------------------------------------------------------------
# 1 · CARGA Y PREPARACIÓN
csv_file = Path("datos/datos_sol_sus_acolite.csv")
df = pd.read_csv(csv_file).drop_duplicates()
df["fecha"] = pd.to_datetime(df["fecha"])

bandas = [c for c in df.columns if c.startswith("B")]
mask   = (df["sol_sus"]>0) & (df[bandas]>0).all(axis=1)
df     = df.loc[mask].copy()

df["log_sol_sus"] = np.log(df["sol_sus"])
for b in bandas: df[f"log_{b}"] = np.log(df[b])

log_cols = [f"log_{b}" for b in bandas]
X_total  = df[log_cols];     y_total = df["log_sol_sus"]

last_date   = df["fecha"].max()
mask_test   = df["fecha"]==last_date
mask_train  = ~mask_test
X_train,y_train = X_total[mask_train],y_total[mask_train]
X_test ,y_test  = X_total[mask_test ],y_total[mask_test ]
grupos_train    = df.loc[mask_train,"fecha"]

# -------------------------------------------------------------
# 2 · AUXILIARES
adj_r2 = lambda r2,n,p: np.nan if n-p-1<=0 else 1-(1-r2)*(n-1)/(n-p-1)
rmse   = lambda a,b: np.sqrt(mean_squared_error(a,b))
def aic_gauss(m,X,y):
    n,k=len(y),X.shape[1]+1
    rss=np.sum((y-m.predict(X))**2)
    return n*np.log(rss/n)+2*k

# -------------------------------------------------------------
# 3 · SELECCIÓN DE BANDA (AIC + LOGO)
cv = LeaveOneGroupOut()
res = []
for b in bandas:
    lb, Xb_tr = f"log_{b}", X_train[[f"log_{b}"]]
    pipe = Pipeline([("s",StandardScaler()),("m",LinearRegression())]).fit(Xb_tr,y_train)
    aic  = aic_gauss(pipe.named_steps["m"],pipe.named_steps["s"].transform(Xb_tr),y_train)
    rm,r2,ra = [],[],[]
    for tr,te in cv.split(Xb_tr,y_train,grupos_train):
        mdl = Pipeline([("s",StandardScaler()),("m",LinearRegression())]).fit(Xb_tr.iloc[tr],y_train.iloc[tr])
        pred = mdl.predict(Xb_tr.iloc[te])
        rm.append(rmse(y_train.iloc[te],pred))
        r2_ = r2_score(y_train.iloc[te],pred)
        r2.append(r2_); ra.append(adj_r2(r2_,len(te),1))
    res.append({"Banda":b,"AIC":aic,
                "RMSE_CV":np.mean(rm),"R2_CV":np.mean(r2),"AdjR2_CV":np.mean(ra)})

tbl = (pd.DataFrame(res).sort_values(["AIC","RMSE_CV"]).round(4).reset_index(drop=True))
display(style_table(tbl,"Ranking por AIC (TRAIN, log)"))

best_band = tbl.loc[0,"Banda"]
display(Markdown(f"**Mejor banda:** `{best_band}`"))

# -------------------------------------------------------------
# 4 · MODELO FINAL
lb_best   = f"log_{best_band}"
pipe_final= Pipeline([("s",StandardScaler()),("m",LinearRegression())]).fit(X_train[[lb_best]],y_train)
a,b_coef  = pipe_final.named_steps["m"].intercept_,pipe_final.named_steps["m"].coef_[0]
display(Markdown("**Ecuación log:**"));  display(Markdown(f"`log(sol_sus) = {a:.4f} + {b_coef:.4f}·log({best_band})`"))
display(Markdown("**Escala real:**"));   display(Markdown(f"`sol_sus = {np.exp(a):.4f} · {best_band}^{b_coef:.4f}`"))

# -------------------------------------------------------------
# 5 · MÉTRICAS TEST (tabla horizontal)
yhat_log  = pipe_final.predict(X_test[[lb_best]])
yhat_lin  = np.exp(yhat_log)
metrics_df = pd.DataFrame({
    "Métrica":["RMSE","R²","R² ajustado"],
    "Log":[rmse(y_test,yhat_log), r2_score(y_test,yhat_log), adj_r2(r2_score(y_test,yhat_log),len(y_test),1)],
    "Real (mg/L)":[rmse(np.exp(y_test),yhat_lin), r2_score(np.exp(y_test),yhat_lin),
                   adj_r2(r2_score(np.exp(y_test),yhat_lin),len(y_test),1)]
})
metrics_df["Log"] = metrics_df["Log"].apply(lambda x:f"{x:.3f}")
metrics_df["Real (mg/L)"]=metrics_df["Real (mg/L)"].apply(lambda x:f"{x:.3f}" if abs(x)<10 else f"{x:.2f}")
display(style_table(metrics_df,"Desempeño en TEST"))

# -------------------------------------------------------------
# 6 · BOOTSTRAP (tabla horizontal)
np.random.seed(123); n_boot=500; fechas=grupos_train.unique()
ib,cb,rm_l,r2_l,ra_l,rm_r,r2_r,ra_r=[],[],[],[],[],[],[],[]
for _ in range(n_boot):
    idx = grupos_train.isin(np.random.choice(fechas,len(fechas),replace=True))
    Xb,yb = X_train[[lb_best]].loc[idx],y_train.loc[idx]
    mdl = Pipeline([("s",StandardScaler()),("m",LinearRegression())]).fit(Xb,yb)
    ib.append(mdl.named_steps["m"].intercept_)
    cb.append(mdl.named_steps["m"].coef_[0])
    pr_log = mdl.predict(X_train[[lb_best]])
    rm_l.append(rmse(y_train,pr_log)); r2_l_ = r2_score(y_train,pr_log); r2_l.append(r2_l_)
    ra_l.append(adj_r2(r2_l_,len(y_train),1))
    pr_lin = np.exp(pr_log)
    rm_r.append(rmse(np.exp(y_train),pr_lin)); r2_r_ = r2_score(np.exp(y_train),pr_lin); r2_r.append(r2_r_)
    ra_r.append(adj_r2(r2_r_,len(y_train),1))

pct = lambda arr: np.percentile(arr,[2.5,97.5])
ci = { "inter":pct(ib), "coef":pct(cb),
       "rm_l":pct(rm_l),"r2_l":pct(r2_l),"ra_l":pct(ra_l),
       "rm_r":pct(rm_r),"r2_r":pct(r2_r),"ra_r":pct(ra_r)}

bs_df = pd.DataFrame({
    "Parámetro / Métrica":[
        "Intercepto",f"Coef {best_band}",
        "RMSE (log)","R² (log)","R² ajustado (log)",
        "RMSE (mg/L)","R² (mg/L)","R² ajustado (mg/L)"
    ],
    "Promedio":[np.mean(ib),np.mean(cb),
                np.mean(rm_l),np.mean(r2_l),np.mean(ra_l),
                np.mean(rm_r),np.mean(r2_r),np.mean(ra_r)],
    "IC 95% inf.":[ci["inter"][0],ci["coef"][0],
                   ci["rm_l"][0],ci["r2_l"][0],ci["ra_l"][0],
                   ci["rm_r"][0],ci["r2_r"][0],ci["ra_r"][0]],
    "IC 95% sup.":[ci["inter"][1],ci["coef"][1],
                   ci["rm_l"][1],ci["r2_l"][1],ci["ra_l"][1],
                   ci["rm_r"][1],ci["r2_r"][1],ci["ra_r"][1]]
})
for col in bs_df.columns[1:]:
    bs_df[col] = bs_df[col].apply(lambda x:f"{x:.4f}" if abs(float(x))<10 else f"{x:.2f}")
display(style_table(bs_df,"Bootstrap (500 réplicas) · IC 95%"))

# -------------------------------------------------------------
# 7 · GRÁFICO Predicho vs Observado  (Train + Test)
#     Izquierda → escala log   ·   Derecha → escala original
plt.rcParams.update({"axes.titlesize":10, "axes.labelsize":8,
                     "xtick.labelsize":7, "ytick.labelsize":7})

TRAIN_CLR = "#9D50A6"   # violeta
TEST_CLR  = "#FFA24B"   # naranja
ID_CLR    = "#17A77E"   # línea identidad

# --- predicciones ---
y_train_log_pred = pipe_final.predict(X_train[[lb_best]])
y_test_log_pred  = pipe_final.predict(X_test [[lb_best]])
y_train_orig, y_train_pred = np.exp(y_train), np.exp(y_train_log_pred)
y_test_orig , y_test_pred  = np.exp(y_test ), np.exp(y_test_log_pred)

fig, ax = plt.subplots(1, 2, figsize=(9, 4), sharex=False, sharey=False)

#  Escala log
ax[0].scatter(y_train, y_train_log_pred, alpha=0.6,
              color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[0].scatter(y_test , y_test_log_pred , alpha=0.6,
              color=TEST_CLR , marker="s", edgecolor="k", label="Test  (■)")

lims = [min(y_train.min(), y_test.min()),
        max(y_train.max(), y_test.max())]
ax[0].plot(lims, lims, "--", color=ID_CLR, lw=2.5, label="y = x")
ax[0].set_xlabel("log(sol_sus) observado")
ax[0].set_ylabel("log(sol_sus) predicho")
ax[0].set_title("Escala log")
ax[0].legend(fontsize=7, frameon=False)
ax[0].grid(ls=":")

eq_log_txt = f"log(sol_sus) = {a:.3f} + {b_coef:.3f}·log({best_band})"
ax[0].text(0.04, 0.96, eq_log_txt, transform=ax[0].transAxes,
           fontsize=7, va="top",
           bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="gray", lw=0.4))

#  Escala original
ax[1].scatter(y_train_orig, y_train_pred, alpha=0.6,
              color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[1].scatter(y_test_orig , y_test_pred , alpha=0.6,
              color=TEST_CLR , marker="s", edgecolor="k", label="Test  (■)")

lims_o = [min(y_train_orig.min(), y_test_orig.min(),
              y_train_pred.min(), y_test_pred.min()),
          max(y_train_orig.max(), y_test_orig.max(),
              y_train_pred.max(), y_test_pred.max())]
ax[1].plot(lims_o, lims_o, "--", color=ID_CLR, lw=2.5, label="y = x")
ax[1].set_xlabel("sol_sus observado (mg/L)")
ax[1].set_ylabel("sol_sus predicho (mg/L)")
ax[1].set_title("Escala original")
ax[1].legend(fontsize=7, frameon=False)
ax[1].grid(ls=":")

eq_orig_txt = f"sol_sus = {np.exp(a):.4f} × {best_band}^{b_coef:.4f}"
ax[1].text(0.04, 0.96, eq_orig_txt, transform=ax[1].transAxes,
           fontsize=7, va="top",
           bbox=dict(boxstyle="round,pad=0.25", fc="white", ec="gray", lw=0.4))

fig.suptitle(f"Predicho vs Observado · Banda {best_band}", fontsize=11, y=1.05)
fig.tight_layout(pad=1.2)
plt.show()

# (opcional) restaurar tamaños por si sigues graficando después
plt.rcParams.update({"axes.titlesize":11, "axes.labelsize":9,
                     "xtick.labelsize":8, "ytick.labelsize":8})

Tabla 7: Ranking por AIC (TRAIN, log)

Banda	AIC	RMSE_CV	R2_CV	AdjR2_CV
B07	-85.218	0.291	-79.607	nan
B08	-84.907	0.285	-77.085	nan
B06	-84.048	0.298	-79.242	nan
B05	-82.491	0.325	-77.426	nan
B8A	-81.287	0.310	-99.788	nan
B04	-79.479	0.340	-96.472	nan
B03	-74.454	0.341	-146.625	nan
B02	-72.724	0.344	-163.983	nan
B01	-70.995	0.342	-163.197	nan
B12	-69.509	0.332	-160.438	nan
B11	-68.415	0.317	-147.291	nan

Mejor banda: B07

Ecuación log:

log(sol_sus) = 4.2330 + 0.1938·log(B07)

Escala real:

sol_sus = 68.9254 · B07^0.1938

Tabla 8: Desempeño en TEST

Métrica	Log	Real (mg/L)
RMSE	0.645	70.77
R²	-1.172	-1.326
R² ajustado	-1.389	-1.558

Tabla 9: Bootstrap (500 réplicas) · IC 95%

Parámetro / Métrica	Promedio	IC 95% inf.	IC 95% sup.
Intercepto	4.2458	4.1789	4.3614
Coef B07	0.1911	0.0600	0.2917
RMSE (log)	0.2440	0.2148	0.3839
R² (log)	0.2398	-0.7610	0.4489
R² ajustado (log)	0.2117	-0.8263	0.4285
RMSE (mg/L)	18.50	15.86	34.59
R² (mg/L)	0.1584	-1.5746	0.4585
R² ajustado (mg/L)	0.1272	-1.6700	0.4384

Conclusión del modelo con validación cruzada

Selecciono la banda B07 como la mejor candidata para estimar los sólidos suspendidos. Esta elección se basa en el valor más bajo del AIC y el menor error promedio en validación cruzada (RMSE_CV). El modelo resultante indica que, a mayor reflectancia en la banda B07, tiende a haber mayor concentración estimada de sólidos suspendidos.

Sin embargo, al evaluar el modelo con los datos más recientes (campaña de prueba), el desempeño cae considerablemente. El R² es negativo, lo que significa que el modelo predice peor que una estimación promedio constante. En otras palabras, el modelo no generaliza bien a nuevos datos.

Para tener una idea más completa de su estabilidad, aplico la técnica de bootstrap con 500 repeticiones. Esto muestra que si bien el coeficiente de la banda B07 es consistentemente positivo (lo que confirma una relación), el rendimiento general del modelo es bajo y muy variable.

Concluyo que, aunque B07 es la mejor banda individual, no alcanza para predecir sólidos suspendidos con precisión. Para mejorar el modelo, es necesario combinar varias bandas, usar índices espectrales o aplicar modelos más complejos.

2.7 Selección de bandas con Forward Selection basado en AIC

En esta sección realizo la selección de un conjunto de bandas espectrales para predecir la concentración de sólidos suspendidos (sol_sus) mediante un modelo de regresión lineal múltiple en escala logarítmica. Utilizo el método de Forward Selection guiado por el criterio de información de Akaike (AIC), que permite incorporar progresivamente aquellas bandas cuya inclusión mejora significativamente el ajuste del modelo sin incrementar de forma innecesaria su complejidad.

En esta sección selecciono un conjunto de bandas espectrales para predecir los sólidos suspendidos (sol_sus) mediante regresión lineal múltiple en escala logarítmica. Utilizo el método de Forward Selection, que comienza con un modelo vacío e incorpora progresivamente las bandas que más reducen el AIC, siempre que la mejora sea significativa (ΔAIC ≥ 1).

Transformo todas las variables a escala logarítmica para linealizar relaciones y reducir la variabilidad. Reservo la última fecha como conjunto de prueba, y entreno el modelo con el resto de los datos. Durante la selección, evalúo el desempeño con validación cruzada agrupada (GroupKFold), que respeta las campañas de muestreo y permite estimar la capacidad del modelo para generalizar a nuevas fechas.

Finalmente, ajusto el modelo con las bandas seleccionadas, calculo métricas en train y test, y realizo un análisis bootstrap para obtener intervalos de confianza de los coeficientes. Esto me permite construir un modelo interpretable, robusto y ajustado a las condiciones reales del monitoreo ambiental.

Código

import numpy as np
import pandas as pd
from pathlib import Path
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import GroupKFold
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
from IPython.display import display, Markdown

# -------------------------------------------------------------
# 1 · CONFIGURACIÓN
DATOS      = Path("datos/datos_sol_sus_acolite.csv")   
N_BOOT     = 2000
DELTA_AIC  = 1.0
RNG_SEED   = 42

# -------------------------------------------------------------
# 2 · ESTILO TABLAS (robusto a versiones pandas)
def style_table(df, caption):
    sty = df.style.format(precision=3)
    if hasattr(sty, "hide"):
        sty = sty.hide(axis="index")
    elif hasattr(sty, "hide_index"):
        sty = sty.hide_index()
    # caption + estilos de alineación
    sty = (sty.set_caption(caption)
              .set_table_styles([
                  {"selector":"caption","props":[("caption-side","top"),
                                                ("font-weight","bold"),
                                                ("margin-bottom","0.5em"),
                                                ("text-align","left")]},
                  {"selector":"th","props":[("text-align","center")]},
                  {"selector":"td","props":[("text-align","right"),
                                            ("padding","0.25em 0.5em")]}
              ], overwrite=False))
    return sty

# -------------------------------------------------------------
# 3 · FUNCIONES AUXILIARES
rmse = lambda a,b: np.sqrt(mean_squared_error(a,b))
aic_g = lambda rss,n,k: n*np.log(rss/n)+2*k
def fit_metrics(Xm, ym):
    m = LinearRegression().fit(Xm, ym)
    pred = m.predict(Xm)
    n,p = Xm.shape
    rss = np.sum((ym - pred)**2)
    return {"AIC":aic_g(rss,n,p+1),
            "RMSE":rmse(ym,pred),
            "R²":r2_score(ym,pred),
            "R²aj":1-(1-r2_score(ym,pred))*(n-1)/(n-p-1),
            "mod":m}

# -------------------------------------------------------------
# 4 · CARGA Y PREPARACIÓN
df = (pd.read_csv(DATOS)
        .sort_values("fecha")
        .reset_index(drop=True))
df["fecha"] = pd.to_datetime(df["fecha"])

bandas = [c for c in df.columns if c.startswith("B")]
mask   = (df["sol_sus"]>0) & (df[bandas]>0).all(axis=1)
df     = df.loc[mask].copy()

df["log_sol_sus"] = np.log(df["sol_sus"])
for b in bandas: df[f"log_{b}"] = np.log(df[b])

X = df[[f"log_{b}" for b in bandas]]; y = df["log_sol_sus"]

last_date = df["fecha"].max()
m_test    = df["fecha"]==last_date
m_train   = ~m_test

X_tr,y_tr = X[m_train],y[m_train]
X_te,y_te = X[m_test ],y[m_test ]
groups_tr = df.loc[m_train,"fecha"]

print(f"TRAIN campañas : {groups_tr.nunique()}")
print(f"TEST  campaña  : {last_date.date()} ({m_test.sum()} muestras)")

# -------------------------------------------------------------
# 5 · FORWARD SELECTION con tabla de candidatos
remaining, selected = list(X_tr.columns), []
rss_null  = np.sum((y_tr - y_tr.mean())**2)
best_aic  = aic_g(rss_null,len(y_tr),1)
hist=[]; paso=0

while remaining:
    cand=[]
    for b in remaining:
        mets = fit_metrics(X_tr[selected+[b]].values, y_tr.values)
        cand.append((mets["AIC"], b, mets))
    cand.sort()

    df_cand = pd.DataFrame({
        "Banda":[b.replace("log_","") for _,b,_ in cand],
        "AIC":[a for a,_,_ in cand]
    })
    df_cand["ΔAIC"] = best_aic - df_cand["AIC"]
    display(Markdown(f"**Paso {paso+1} - candidatos:**"))
    display(style_table(df_cand,""))

    new_aic,best_band,best_met = cand[0]
    if best_aic - new_aic >= DELTA_AIC:
        selected.append(best_band); remaining.remove(best_band)
        best_aic=new_aic; paso+=1
        hist.append({"Paso":paso,
                     "Bandas":", ".join(c.replace("log_","") for c in selected),
                     "AIC":best_aic,
                     "RMSE":best_met["RMSE"],
                     "R²":best_met["R²"],
                     "R²aj":best_met["R²aj"]})
    else:
        break

if not selected:
    selected=[cand[0][1]]; best_aic=cand[0][0]; best_met=cand[0][2]
    hist=[{"Paso":1,"Bandas":selected[0].replace("log_",""),
           "AIC":best_aic,"RMSE":best_met["RMSE"],
           "R²":best_met["R²"],"R²aj":best_met["R²aj"]}]

hist_df=pd.DataFrame(hist)[["Paso","Bandas","AIC","RMSE","R²","R²aj"]]
display(style_table(hist_df,"Historial forward‑AIC"))

sel_log=selected; sel=[c.replace("log_","") for c in sel_log]
display(Markdown(f"**Bandas finales:** {', '.join(sel)}"))

# -------------------------------------------------------------
# 6 · GroupKFold CV
gkf=GroupKFold(n_splits=min(5,groups_tr.nunique()))
rmse_cv=[rmse(y_tr.iloc[te],
              LinearRegression().fit(X_tr.iloc[tr][sel_log],y_tr.iloc[tr])
              .predict(X_tr.iloc[te][sel_log]))
         for tr,te in gkf.split(X_tr[sel_log],y_tr,groups_tr)]
display(Markdown(f"**GroupKFold CV (log) RMSE:** {np.mean(rmse_cv):.3f} ± {np.std(rmse_cv):.3f}"))

# -------------------------------------------------------------
# 7 · MODELO FINAL + MÉTRICAS
mod = LinearRegression().fit(X_tr[sel_log], y_tr)
pred_tr_log = mod.predict(X_tr[sel_log])
pred_te_log = mod.predict(X_te[sel_log])

y_tr_lin, y_te_lin = np.exp(y_tr), np.exp(y_te)
pred_tr_lin, pred_te_lin = np.exp(pred_tr_log), np.exp(pred_te_log)

metrics_df = pd.DataFrame({
    "Métrica": ["RMSE", "R²"],
    "Train log": [rmse(y_tr, pred_tr_log), r2_score(y_tr, pred_tr_log)],
    "Test log": [rmse(y_te, pred_te_log), r2_score(y_te, pred_te_log)],
    "Train orig": [rmse(y_tr_lin, pred_tr_lin), r2_score(y_tr_lin, pred_tr_lin)],
    "Test orig": [rmse(y_te_lin, pred_te_lin), r2_score(y_te_lin, pred_te_lin)]
})
for col in metrics_df.columns[1:]:
    metrics_df[col] = metrics_df[col].apply(lambda x: f"{x:.3f}")

display(style_table(metrics_df, "Desempeño en Train y Test"))

# -------------------------------------------------------------
# 8 · BOOTSTRAP coeficientes
rng = np.random.default_rng(RNG_SEED)
dates = groups_tr.unique()
coef_bs = []

for _ in range(N_BOOT):
    idx = groups_tr.isin(rng.choice(dates, len(dates), replace=True))
    coef_bs.append(LinearRegression().fit(X_tr.loc[idx, sel_log], y_tr.loc[idx]).coef_)

coef_bs = np.array(coef_bs)
mean = coef_bs.mean(axis=0)
lo, hi = np.percentile(coef_bs, [2.5, 97.5], axis=0)

bs_df = pd.DataFrame({"Banda": sel, "Media β": mean, "IC inf": lo, "IC sup": hi})
display(style_table(bs_df, "Bootstrap (coeficientes, 95 % CI)"))

# -------------------------------------------------------------
# 9 · ECUACIONES Y BANDAS SELECCIONADAS
alpha = np.exp(mod.intercept_)
expr_log = " + ".join(f"{c:.4f}·log({b})" for c, b in zip(mod.coef_, sel))
expr_orig = " · ".join(f"{b}^{c:.4f}" for c, b in zip(mod.coef_, sel))

display(Markdown("**Bandas seleccionadas:**"))
display(Markdown(f"`{', '.join(sel)}`"))

display(Markdown("**Ecuación en escala logarítmica:**"))
display(Markdown(f"`log(sol_sus) = {mod.intercept_:.4f} + {expr_log}`"))

display(Markdown("**Ecuación en escala original:**"))
display(Markdown(f"`sol_sus = {alpha:.4f} · {expr_orig}`"))

# -------------------------------------------------------------
# 10 · GRÁFICO Predicho vs Observado (Train+Test)
plt.rcParams.update({"axes.titlesize":10,"axes.labelsize":8,
                     "xtick.labelsize":7,"ytick.labelsize":7})

TRAIN_CLR = "#9D50A6"
TEST_CLR = "#FFA24B"
ID_CLR = "#17A77E"

fig, ax = plt.subplots(1, 2, figsize=(9, 4))
# Escala log
ax[0].scatter(y_tr, pred_tr_log, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[0].scatter(y_te, pred_te_log, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims = [min(y_tr.min(), y_te.min()), max(y_tr.max(), y_te.max())]
ax[0].plot(lims, lims, "--", color=ID_CLR, lw=2.5, label="y = x")
ax[0].set(xlabel="log(sol_sus) obs.", ylabel="log(sol_sus) pred.", title="Escala log")
ax[0].legend(fontsize=7, frameon=False)
ax[0].grid(ls=":")

# Escala original
ax[1].scatter(y_tr_lin, pred_tr_lin, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[1].scatter(y_te_lin, pred_te_lin, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims_o = [min(y_tr_lin.min(), y_te_lin.min(), pred_tr_lin.min(), pred_te_lin.min()),
          max(y_tr_lin.max(), y_te_lin.max(), pred_tr_lin.max(), pred_te_lin.max())]
ax[1].plot(lims_o, lims_o, "--", color=ID_CLR, lw=2.5)
ax[1].set(xlabel="sol_sus obs. (mg/L)", ylabel="sol_sus pred. (mg/L)", title="Escala original")
ax[1].legend(fontsize=7, frameon=False)
ax[1].grid(ls=":")

fig.suptitle(f"Predicho vs Observado · ΔAIC ≥ {DELTA_AIC}", fontsize=11, y=1.05)
fig.tight_layout(pad=1.2)
plt.show()

# restaurar estilos si continuás graficando
plt.rcParams.update({"axes.titlesize":11,"axes.labelsize":9,
                     "xtick.labelsize":8,"ytick.labelsize":8})

TRAIN campañas : 5
TEST  campaña  : 2025-06-09 (12 muestras)

Paso 1 - candidatos:

Banda	AIC	ΔAIC
B07	-85.218	15.280
B08	-84.907	14.970
B06	-84.048	14.110
B05	-82.491	12.553
B8A	-81.287	11.350
B04	-79.479	9.542
B03	-74.454	4.517
B02	-72.724	2.787
B01	-70.994	1.057
B12	-69.509	-0.429
B11	-68.415	-1.523

Paso 2 - candidatos:

Banda	AIC	ΔAIC
B02	-108.017	22.799
B01	-99.457	14.240
B03	-97.017	11.799
B11	-95.445	10.228
B12	-88.549	3.332
B8A	-87.786	2.568
B06	-85.179	-0.039
B04	-84.967	-0.251
B05	-83.319	-1.899
B08	-83.234	-1.984

Paso 3 - candidatos:

Banda	AIC	ΔAIC
B11	-108.355	0.338
B08	-108.291	0.274
B03	-107.375	-0.641
B8A	-107.364	-0.653
B04	-107.084	-0.933
B05	-106.992	-1.025
B12	-106.731	-1.286
B01	-106.713	-1.304
B06	-106.021	-1.996

Tabla 10: Historial forward‑AIC

Paso	Bandas	AIC	RMSE	R²	R²aj
1	B07	-85.218	0.215	0.449	0.429
2	B07, B02	-108.017	0.140	0.766	0.748

Bandas finales: B07, B02

GroupKFold CV (log) RMSE: 0.125 ± 0.058

Tabla 11: Desempeño en Train y Test

Métrica	Train log	Test log	Train orig	Test orig
RMSE	0.140	0.314	9.209	36.331
R²	0.766	0.484	0.817	0.387

Tabla 12: Bootstrap (coeficientes, 95 % CI)

Banda	Media β	IC inf	IC sup
B07	0.635	0.429	0.713
B02	-0.698	-0.881	-0.141

Bandas seleccionadas:

B07, B02

Ecuación en escala logarítmica:

log(sol_sus) = 3.9813 + 0.6438·log(B07) + -0.7532·log(B02)

Ecuación en escala original:

sol_sus = 53.5891 · B07^0.6438 · B02^-0.7532

Conclusión del modelo seleccionado

Selecciono las bandas B07 y B02 como las mejores para estimar los sólidos suspendidos, basándome en la reducción significativa del AIC y la mejora en el error promedio de validación cruzada (RMSE_CV). El modelo indica que la reflectancia en B07 tiene un efecto positivo en la concentración estimada, mientras que B02 aporta un efecto negativo, lo que puede relacionarse con propiedades ópticas del agua y sedimentos.

Al evaluar el modelo con datos nuevos (campaña de prueba), el desempeño es razonablemente bueno, con un R² positivo que refleja una capacidad aceptable de generalización. El análisis bootstrap confirma la estabilidad de los coeficientes y las métricas, entregando confianza en la robustez del modelo.

En resumen, combinar estas dos bandas mejora sustancialmente la predicción en comparación con usar una sola banda, pero para lograr mayor precisión probablemente sea necesario explorar índices espectrales o modelos más avanzados.

2.8 Selección de variables mediante Forward Selection con AIC usando razones simples de bandas

En esta etapa exploro combinaciones de variables formadas por razones simples entre bandas espectrales, todas en escala logarítmica, para mejorar la predicción de los sólidos suspendidos (sol_sus). Aplico un método de Forward Selection guiado por el criterio de información de Akaike (AIC), que me permite construir un modelo parsimonioso incorporando solo aquellas razones que aportan mejoras significativas en el ajuste.

El análisis se realiza sobre el conjunto de entrenamiento definido previamente, manteniendo la partición temporal para test y respetando la estructura grupal en la validación cruzada. Además, implemento un análisis bootstrap para estimar la incertidumbre en los coeficientes y evaluar la estabilidad del modelo.

Con este procedimiento busco identificar relaciones espectrales más complejas que mejor expliquen la variabilidad de sólidos suspendidos, maximizando la capacidad predictiva sin sobreajustar, para un monitoreo ambiental más preciso y confiable.

Código

# 0 · LIBRERÍAS
import numpy as np
import pandas as pd
from pathlib import Path
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import GroupKFold
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
from IPython.display import display, Markdown

# -------------------------------------------------------------
# 1 · CONFIGURACIÓN
DATOS     = Path("datos/datos_sol_sus_acolite.csv")     # ← ajustá si hace falta
DELTA_AIC = 2                                            # ΔAIC mínimo
N_BOOT    = 2000
RNG_SEED  = 42

# -------------------------------------------------------------
# 2 · FUNCIONES AUXILIARES
def aic_gauss(rss, n, k):           # k = nº coef + intercepto
    return n * np.log(rss / n) + 2 * k

def fit_metrics(X, y):
    """Ajusta LR y devuelve métricas + modelo."""
    m    = LinearRegression().fit(X, y)
    pred = m.predict(X)
    rss  = np.sum((y - pred) ** 2)
    n, p = X.shape
    return {
        "AIC" : aic_gauss(rss, n, p + 1),
        "RMSE": np.sqrt(mean_squared_error(y, pred)),
        "R²"  : r2_score(y, pred),
        "R²aj": 1 - (1 - r2_score(y, pred)) * (n - 1) / (n - p - 1),
        "mod" : m
    }

def rmse(a, b): return np.sqrt(mean_squared_error(a, b))

# -------------------------------------------------------------
# 3 · CARGA Y PREPARACIÓN DE DATOS
df = pd.read_csv(DATOS)
df["fecha"] = pd.to_datetime(df["fecha"])

bandas = [c for c in df.columns if c.startswith("B")]
mask   = (df["sol_sus"] > 0) & (df[bandas] > 0).all(axis=1)
df     = df.loc[mask].copy()

# variable respuesta en log
df["log_sol_sus"] = np.log(df["sol_sus"])

# -------------------------------------------------------------
# 3.1 · GENERAR RAZONES SIMPLES DE BANDAS EN LOGARITMO
ratio_cols = []
for i, b1 in enumerate(bandas):
    for b2 in bandas[i+1:]:
        col = f"log({b1}/{b2})"
        df[col] = np.log(df[b1]) - np.log(df[b2])
        ratio_cols.append(col)

X = df[ratio_cols]                # predictores = razones simples en log
y = df["log_sol_sus"]

# -------------------------------------------------------------
# 3.2 · PARTICIÓN TEMPORAL
last_date = df["fecha"].max()
m_test    = df["fecha"] == last_date
m_train   = ~m_test

X_tr, y_tr = X.loc[m_train], y.loc[m_train]
X_te, y_te = X.loc[m_test] , y.loc[m_test]
groups_tr  = df.loc[m_train, "fecha"]

print(f"TRAIN campañas : {groups_tr.nunique()}")
print(f"TEST  campaña  : {last_date.date()} ({m_test.sum()} muestras)")

# -------------------------------------------------------------
# 4 · FORWARD SELECTION (razones simples logarítmicas) con AIC + Δ
remaining = ratio_cols.copy()
selected  = []

# modelo nulo
rss_null  = np.sum((y_tr - y_tr.mean())**2)
best_aic  = aic_gauss(rss_null, len(y_tr), k=1)

hist = []
paso = 0
while remaining:
    candidatos = []
    for col in remaining:
        cols = selected + [col]
        m    = fit_metrics(X_tr[cols].values, y_tr.values)
        candidatos.append((m["AIC"], col, m))
    candidatos.sort()

    new_aic, best_col, best_met = candidatos[0]
    if best_aic - new_aic >= DELTA_AIC:     # mejora ≥ ΔAIC
        selected.append(best_col)
        remaining.remove(best_col)
        best_aic = new_aic
        paso += 1
        hist.append({
            "Paso"  : paso,
            "Bandas": ", ".join(selected),
            "AIC"   : best_aic,
            "RMSE"  : best_met["RMSE"],
            "R²"    : best_met["R²"],
            "R²aj"  : best_met["R²aj"]
        })
    else:
        break

# fallback (si nada mejora)
if not selected:
    selected = [min(candidatos, key=lambda t: t[0])[1]]
    best_aic = min(candidatos)[0]
    best_met = min(candidatos)[2]
    hist.append({
        "Paso"  : 1,
        "Bandas": selected[0],
        "AIC"   : best_aic,
        "RMSE"  : best_met["RMSE"],
        "R²"    : best_met["R²"],
        "R²aj"  : best_met["R²aj"]
    })

hist_df = (pd.DataFrame(hist)
             .loc[:, ["Paso", "Bandas", "AIC", "RMSE", "R²", "R²aj"]])

# Mostrar tabla historial con caption y sin índice
display(
    hist_df.style
    .format({"AIC": "{:.3f}", "RMSE": "{:.3f}", "R²": "{:.3f}", "R²aj": "{:.3f}"})
    .set_caption("Historial de selección")
    .hide(axis="index")
)

# Mostrar bandas seleccionadas en negrita Markdown
display(Markdown(f"**Bandas seleccionadas:**  `{', '.join(selected)}`"))

# -------------------------------------------------------------
# 5 · VALIDACIÓN GROUP‑KFOLD
gkf = GroupKFold(n_splits=min(5, groups_tr.nunique()))
rmse_cv = []
for tr_idx, te_idx in gkf.split(X_tr[selected], y_tr, groups_tr):
    m    = LinearRegression().fit(X_tr.iloc[tr_idx][selected], y_tr.iloc[tr_idx])
    pred = m.predict(X_tr.iloc[te_idx][selected])
    rmse_cv.append(rmse(y_tr.iloc[te_idx], pred))

# Mostrar GroupKFold CV en negrita Markdown
display(Markdown(f"**GroupKFold CV (log) RMSE:** {np.mean(rmse_cv):.3f} ± {np.std(rmse_cv):.3f}"))

# -------------------------------------------------------------
# 6 · MODELO FINAL Y MÉTRICAS
mod = LinearRegression().fit(X_tr[selected], y_tr)

# TRAIN
pred_tr_log = mod.predict(X_tr[selected])
rmse_tr_log = rmse(y_tr, pred_tr_log)
r2_tr_log   = r2_score(y_tr, pred_tr_log)

# TEST
pred_te_log = mod.predict(X_te[selected])
rmse_te_log = rmse(y_te, pred_te_log)
r2_te_log   = r2_score(y_te, pred_te_log)

# escala original
y_tr_lin, y_te_lin = np.exp(y_tr), np.exp(y_te)
pred_tr_lin, pred_te_lin = np.exp(pred_tr_log), np.exp(pred_te_log)

rmse_tr_lin = rmse(y_tr_lin, pred_tr_lin)
r2_tr_lin   = r2_score(y_tr_lin, pred_tr_lin)
rmse_te_lin = rmse(y_te_lin, pred_te_lin)
r2_te_lin   = r2_score(y_te_lin, pred_te_lin)

metrics_df = pd.DataFrame({
    "Métrica"    : ["RMSE", "R²"],
    "Train log"  : [rmse_tr_log, r2_tr_log],
    "Test log"   : [rmse_te_log, r2_te_log],
    "Train orig" : [rmse_tr_lin, r2_tr_lin],
    "Test orig"  : [rmse_te_lin, r2_te_lin]
})

# Mostrar tabla métricas con caption y sin índice
display(
    metrics_df.style
    .format({
        "Train log": "{:.3f}", "Test log": "{:.3f}",
        "Train orig": "{:.3f}", "Test orig": "{:.3f}"
    })
    .set_caption("Métricas finales")
    .hide(axis="index")
)

# -------------------------------------------------------------
# 7 · BOOTSTRAP AGRUPADO (coef 95 % CI)
rng    = np.random.default_rng(RNG_SEED)
dates  = groups_tr.unique()
coef_bs = []

for _ in range(N_BOOT):
    samp = rng.choice(dates, len(dates), replace=True)
    idx  = groups_tr.isin(samp)
    m_bs = LinearRegression().fit(X_tr.loc[idx, selected], y_tr.loc[idx])
    coef_bs.append(m_bs.coef_)

coef_bs = np.array(coef_bs)
mean_bs = coef_bs.mean(axis=0)
ci_bs   = np.percentile(coef_bs, [2.5, 97.5], axis=0)

bs_df = pd.DataFrame({
    "Banda"       : selected,
    "Media β"     : mean_bs,
    "IC 95% inf." : ci_bs[0],
    "IC 95% sup." : ci_bs[1]
})

# Mostrar tabla bootstrap con caption y sin índice
display(
    bs_df.style
    .format({
        "Media β": "{:.4f}",
        "IC 95% inf.": "{:.4f}",
        "IC 95% sup.": "{:.4f}"
    })
    .set_caption("Bootstrap (coeficientes, 95 % CI)")
    .hide(axis="index")
)

# -------------------------------------------------------------
# 8 · ECUACIONES
alpha      = np.exp(mod.intercept_)
expr_log   = " + ".join([f"{c:.4f}·{r}" for c, r in zip(mod.coef_, selected)])
expr_orig  = " · ".join([f"{r}^{c:.4f}" for c, r in zip(mod.coef_, selected)])

display(Markdown("Ecuación (log)"))
display(Markdown(f"`log(sol_sus) = {mod.intercept_:.4f} + {expr_log}`"))

display(Markdown("Ecuación (original)"))
display(Markdown(f"`sol_sus = {alpha:.4f} · {expr_orig}`"))

# -------------------------------------------------------------
# 9 · GRÁFICAS DE PARIDAD (TRAIN + TEST) - estilo especificado
plt.rcParams.update({"axes.titlesize":10,"axes.labelsize":8,
                     "xtick.labelsize":7,"ytick.labelsize":7})

TRAIN_CLR = "#9D50A6"
TEST_CLR = "#FFA24B"
ID_CLR = "#17A77E"

fig, ax = plt.subplots(1, 2, figsize=(9, 4))
# Escala log
ax[0].scatter(y_tr, pred_tr_log, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[0].scatter(y_te, pred_te_log, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims = [min(y_tr.min(), y_te.min()), max(y_tr.max(), y_te.max())]
ax[0].plot(lims, lims, "--", color=ID_CLR, lw=2.5, label="y = x")
ax[0].set(xlabel="log(sol_sus) obs.", ylabel="log(sol_sus) pred.", title="Escala log")
ax[0].legend(fontsize=7, frameon=False)
ax[0].grid(ls=":")

# Escala original
ax[1].scatter(y_tr_lin, pred_tr_lin, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[1].scatter(y_te_lin, pred_te_lin, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims_o = [min(y_tr_lin.min(), y_te_lin.min(), pred_tr_lin.min(), pred_te_lin.min()),
          max(y_tr_lin.max(), y_te_lin.max(), pred_tr_lin.max(), pred_te_lin.max())]
ax[1].plot(lims_o, lims_o, "--", color=ID_CLR, lw=2.5)
ax[1].set(xlabel="sol_sus obs. (mg/L)", ylabel="sol_sus pred. (mg/L)", title="Escala original")
ax[1].legend(fontsize=7, frameon=False)
ax[1].grid(ls=":")

fig.suptitle(f"Predicho vs Observado · ΔAIC ≥ {DELTA_AIC}", fontsize=11, y=1.05)
fig.tight_layout(pad=1.2)
plt.show()

# restaurar estilos si continuás graficando
plt.rcParams.update({"axes.titlesize":11,"axes.labelsize":9,
                     "xtick.labelsize":8,"ytick.labelsize":8})

TRAIN campañas : 5
TEST  campaña  : 2025-06-09 (12 muestras)

Tabla 13: Historial de selección

Paso	Bandas	AIC	RMSE	R²	R²aj
1	log(B02/B07)	-106.889	0.148	0.739	0.729
2	log(B02/B07), log(B05/B11)	-110.321	0.135	0.784	0.767

Bandas seleccionadas: log(B02/B07), log(B05/B11)

GroupKFold CV (log) RMSE: 0.119 ± 0.048

Tabla 14: Métricas finales

Métrica	Train log	Test log	Train orig	Test orig
RMSE	0.135	0.388	9.258	45.927
R²	0.784	0.212	0.816	0.021

Tabla 15: Bootstrap (coeficientes, 95 % CI)

Banda	Media β	IC 95% inf.	IC 95% sup.
log(B02/B07)	-0.5842	-0.6164	-0.5265
log(B05/B11)	0.0239	-0.0012	0.0348

Ecuación (log)

log(sol_sus) = 4.2142 + -0.5853·log(B02/B07) + 0.0256·log(B05/B11)

Ecuación (original)

sol_sus = 67.6397 · log(B02/B07)^-0.5853 · log(B05/B11)^0.0256

Conclusión del modelo seleccionado

En esta etapa, al usar razones simples de bandas espectrales en lugar de bandas individuales, logro un modelo con mejor ajuste en entrenamiento, con un R² de 0.784 frente a 0.766 del modelo previo. También el error promedio en la validación cruzada mejora ligeramente (RMSE CV 0.119 vs. 0.125), lo que indica una selección más eficiente de variables.

Sin embargo, al evaluar con los datos de la campaña más reciente, el desempeño baja significativamente: el R² cae a 0.212, mucho menor que el 0.484 obtenido anteriormente, y el RMSE aumenta, mostrando que la capacidad de generalización sigue siendo limitada.

El análisis bootstrap confirma que la razón log(B02/B07) tiene un coeficiente negativo estable y significativo, mientras que log(B05/B11) aporta muy poco y su coeficiente incluye cero en el intervalo de confianza, indicando que su efecto es marginal.

Desde el punto de vista físico, el modelo indica que a medida que aumenta la razón B02/B07, la concentración de sólidos suspendidos tiende a disminuir, mientras que la razón B05/B11 tiene un efecto muy débil e incierto.

A pesar de que el uso de razones simples resulta en un modelo más compacto y con una mejora modesta en el ajuste, la capacidad para predecir sólidos suspendidos en campañas futuras continúa siendo limitada.

2.8.1 Ajuste del umbral ΔAIC para selección de variables

Para evaluar la robustez y simplicidad del modelo, aumento el umbral mínimo de mejora en el criterio AIC (ΔAIC) de 2 a 4. Esto implica que la selección hacia adelante solo incorpora nuevas variables si la reducción en AIC es al menos 4, promoviendo modelos más parsimoniosos y evitando incluir variables con mejoras marginales que podrían no generalizar bien.

Este ajuste busca balancear mejor el compromiso entre ajuste y complejidad, reduciendo el riesgo de sobreajuste, especialmente relevante al aplicar el modelo a campañas de muestreo futuras.

Código

# 0 · LIBRERÍAS
import numpy as np
import pandas as pd
from pathlib import Path
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import GroupKFold
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
from IPython.display import display, Markdown

# -------------------------------------------------------------
# 1 · CONFIGURACIÓN
DATOS     = Path("datos/datos_sol_sus_acolite.csv")     # ← ajustá si hace falta
DELTA_AIC = 4                                            # ΔAIC mínimo
N_BOOT    = 2000
RNG_SEED  = 42

# -------------------------------------------------------------
# 2 · FUNCIONES AUXILIARES
def aic_gauss(rss, n, k):           # k = nº coef + intercepto
    return n * np.log(rss / n) + 2 * k

def fit_metrics(X, y):
    """Ajusta LR y devuelve métricas + modelo."""
    m    = LinearRegression().fit(X, y)
    pred = m.predict(X)
    rss  = np.sum((y - pred) ** 2)
    n, p = X.shape
    return {
        "AIC" : aic_gauss(rss, n, p + 1),
        "RMSE": np.sqrt(mean_squared_error(y, pred)),
        "R²"  : r2_score(y, pred),
        "R²aj": 1 - (1 - r2_score(y, pred)) * (n - 1) / (n - p - 1),
        "mod" : m
    }

def rmse(a, b): return np.sqrt(mean_squared_error(a, b))

# -------------------------------------------------------------
# 3 · CARGA Y PREPARACIÓN DE DATOS
df = pd.read_csv(DATOS)
df["fecha"] = pd.to_datetime(df["fecha"])

bandas = [c for c in df.columns if c.startswith("B")]
mask   = (df["sol_sus"] > 0) & (df[bandas] > 0).all(axis=1)
df     = df.loc[mask].copy()

# variable respuesta en log
df["log_sol_sus"] = np.log(df["sol_sus"])

# -------------------------------------------------------------
# 3.1 · GENERAR RAZONES SIMPLES DE BANDAS EN LOGARITMO
ratio_cols = []
for i, b1 in enumerate(bandas):
    for b2 in bandas[i+1:]:
        col = f"log({b1}/{b2})"
        df[col] = np.log(df[b1]) - np.log(df[b2])
        ratio_cols.append(col)

X = df[ratio_cols]                # predictores = razones simples en log
y = df["log_sol_sus"]

# -------------------------------------------------------------
# 3.2 · PARTICIÓN TEMPORAL
last_date = df["fecha"].max()
m_test    = df["fecha"] == last_date
m_train   = ~m_test

X_tr, y_tr = X.loc[m_train], y.loc[m_train]
X_te, y_te = X.loc[m_test] , y.loc[m_test]
groups_tr  = df.loc[m_train, "fecha"]

print(f"TRAIN campañas : {groups_tr.nunique()}")
print(f"TEST  campaña  : {last_date.date()} ({m_test.sum()} muestras)")

# -------------------------------------------------------------
# 4 · FORWARD SELECTION (razones simples logarítmicas) con AIC + Δ
remaining = ratio_cols.copy()
selected  = []

# modelo nulo
rss_null  = np.sum((y_tr - y_tr.mean())**2)
best_aic  = aic_gauss(rss_null, len(y_tr), k=1)

hist = []
paso = 0
while remaining:
    candidatos = []
    for col in remaining:
        cols = selected + [col]
        m    = fit_metrics(X_tr[cols].values, y_tr.values)
        candidatos.append((m["AIC"], col, m))
    candidatos.sort()

    new_aic, best_col, best_met = candidatos[0]
    if best_aic - new_aic >= DELTA_AIC:     # mejora ≥ ΔAIC
        selected.append(best_col)
        remaining.remove(best_col)
        best_aic = new_aic
        paso += 1
        hist.append({
            "Paso"  : paso,
            "Bandas": ", ".join(selected),
            "AIC"   : best_aic,
            "RMSE"  : best_met["RMSE"],
            "R²"    : best_met["R²"],
            "R²aj"  : best_met["R²aj"]
        })
    else:
        break

# fallback (si nada mejora)
if not selected:
    selected = [min(candidatos, key=lambda t: t[0])[1]]
    best_aic = min(candidatos)[0]
    best_met = min(candidatos)[2]
    hist.append({
        "Paso"  : 1,
        "Bandas": selected[0],
        "AIC"   : best_aic,
        "RMSE"  : best_met["RMSE"],
        "R²"    : best_met["R²"],
        "R²aj"  : best_met["R²aj"]
    })

hist_df = (pd.DataFrame(hist)
             .loc[:, ["Paso", "Bandas", "AIC", "RMSE", "R²", "R²aj"]])

# Mostrar tabla historial con caption y sin índice
display(
    hist_df.style
    .format({"AIC": "{:.3f}", "RMSE": "{:.3f}", "R²": "{:.3f}", "R²aj": "{:.3f}"})
    .set_caption("Historial de selección")
    .hide(axis="index")
)

# Mostrar bandas seleccionadas en negrita Markdown
display(Markdown(f"**Bandas seleccionadas:**  `{', '.join(selected)}`"))

# -------------------------------------------------------------
# 5 · VALIDACIÓN GROUP‑KFOLD
gkf = GroupKFold(n_splits=min(5, groups_tr.nunique()))
rmse_cv = []
for tr_idx, te_idx in gkf.split(X_tr[selected], y_tr, groups_tr):
    m    = LinearRegression().fit(X_tr.iloc[tr_idx][selected], y_tr.iloc[tr_idx])
    pred = m.predict(X_tr.iloc[te_idx][selected])
    rmse_cv.append(rmse(y_tr.iloc[te_idx], pred))

# Mostrar GroupKFold CV en negrita Markdown
display(Markdown(f"**GroupKFold CV (log) RMSE:** {np.mean(rmse_cv):.3f} ± {np.std(rmse_cv):.3f}"))

# -------------------------------------------------------------
# 6 · MODELO FINAL Y MÉTRICAS
mod = LinearRegression().fit(X_tr[selected], y_tr)

# TRAIN
pred_tr_log = mod.predict(X_tr[selected])
rmse_tr_log = rmse(y_tr, pred_tr_log)
r2_tr_log   = r2_score(y_tr, pred_tr_log)

# TEST
pred_te_log = mod.predict(X_te[selected])
rmse_te_log = rmse(y_te, pred_te_log)
r2_te_log   = r2_score(y_te, pred_te_log)

# escala original
y_tr_lin, y_te_lin = np.exp(y_tr), np.exp(y_te)
pred_tr_lin, pred_te_lin = np.exp(pred_tr_log), np.exp(pred_te_log)

rmse_tr_lin = rmse(y_tr_lin, pred_tr_lin)
r2_tr_lin   = r2_score(y_tr_lin, pred_tr_lin)
rmse_te_lin = rmse(y_te_lin, pred_te_lin)
r2_te_lin   = r2_score(y_te_lin, pred_te_lin)

metrics_df = pd.DataFrame({
    "Métrica"    : ["RMSE", "R²"],
    "Train log"  : [rmse_tr_log, r2_tr_log],
    "Test log"   : [rmse_te_log, r2_te_log],
    "Train orig" : [rmse_tr_lin, r2_tr_lin],
    "Test orig"  : [rmse_te_lin, r2_te_lin]
})

# Mostrar tabla métricas con caption y sin índice
display(
    metrics_df.style
    .format({
        "Train log": "{:.3f}", "Test log": "{:.3f}",
        "Train orig": "{:.3f}", "Test orig": "{:.3f}"
    })
    .set_caption("Métricas finales")
    .hide(axis="index")
)

# -------------------------------------------------------------
# 7 · BOOTSTRAP AGRUPADO (coef 95 % CI)
rng    = np.random.default_rng(RNG_SEED)
dates  = groups_tr.unique()
coef_bs = []

for _ in range(N_BOOT):
    samp = rng.choice(dates, len(dates), replace=True)
    idx  = groups_tr.isin(samp)
    m_bs = LinearRegression().fit(X_tr.loc[idx, selected], y_tr.loc[idx])
    coef_bs.append(m_bs.coef_)

coef_bs = np.array(coef_bs)
mean_bs = coef_bs.mean(axis=0)
ci_bs   = np.percentile(coef_bs, [2.5, 97.5], axis=0)

bs_df = pd.DataFrame({
    "Banda"       : selected,
    "Media β"     : mean_bs,
    "IC 95% inf." : ci_bs[0],
    "IC 95% sup." : ci_bs[1]
})

# Mostrar tabla bootstrap con caption y sin índice
display(
    bs_df.style
    .format({
        "Media β": "{:.4f}",
        "IC 95% inf.": "{:.4f}",
        "IC 95% sup.": "{:.4f}"
    })
    .set_caption("Bootstrap (coeficientes, 95 % CI)")
    .hide(axis="index")
)

# -------------------------------------------------------------
# 8 · ECUACIONES
alpha      = np.exp(mod.intercept_)
expr_log   = " + ".join([f"{c:.4f}·{r}" for c, r in zip(mod.coef_, selected)])
expr_orig  = " · ".join([f"{r}^{c:.4f}" for c, r in zip(mod.coef_, selected)])

display(Markdown("Ecuación (log)"))
display(Markdown(f"`log(sol_sus) = {mod.intercept_:.4f} + {expr_log}`"))

display(Markdown("Ecuación (original)"))
display(Markdown(f"`sol_sus = {alpha:.4f} · {expr_orig}`"))

# -------------------------------------------------------------
# 9 · GRÁFICAS DE PARIDAD (TRAIN + TEST) - estilo especificado
plt.rcParams.update({"axes.titlesize":10,"axes.labelsize":8,
                     "xtick.labelsize":7,"ytick.labelsize":7})

TRAIN_CLR = "#9D50A6"
TEST_CLR = "#FFA24B"
ID_CLR = "#17A77E"

fig, ax = plt.subplots(1, 2, figsize=(9, 4))
# Escala log
ax[0].scatter(y_tr, pred_tr_log, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[0].scatter(y_te, pred_te_log, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims = [min(y_tr.min(), y_te.min()), max(y_tr.max(), y_te.max())]
ax[0].plot(lims, lims, "--", color=ID_CLR, lw=2.5, label="y = x")
ax[0].set(xlabel="log(sol_sus) obs.", ylabel="log(sol_sus) pred.", title="Escala log")
ax[0].legend(fontsize=7, frameon=False)
ax[0].grid(ls=":")

# Escala original
ax[1].scatter(y_tr_lin, pred_tr_lin, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[1].scatter(y_te_lin, pred_te_lin, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims_o = [min(y_tr_lin.min(), y_te_lin.min(), pred_tr_lin.min(), pred_te_lin.min()),
          max(y_tr_lin.max(), y_te_lin.max(), pred_tr_lin.max(), pred_te_lin.max())]
ax[1].plot(lims_o, lims_o, "--", color=ID_CLR, lw=2.5)
ax[1].set(xlabel="sol_sus obs. (mg/L)", ylabel="sol_sus pred. (mg/L)", title="Escala original")
ax[1].legend(fontsize=7, frameon=False)
ax[1].grid(ls=":")

fig.suptitle(f"Predicho vs Observado · ΔAIC ≥ {DELTA_AIC}", fontsize=11, y=1.05)
fig.tight_layout(pad=1.2)
plt.show()

# restaurar estilos si continuás graficando
plt.rcParams.update({"axes.titlesize":11,"axes.labelsize":9,
                     "xtick.labelsize":8,"ytick.labelsize":8})

TRAIN campañas : 5
TEST  campaña  : 2025-06-09 (12 muestras)

Tabla 16: Historial de selección

Paso	Bandas	AIC	RMSE	R²	R²aj
1	log(B02/B07)	-106.889	0.148	0.739	0.729

Bandas seleccionadas: log(B02/B07)

GroupKFold CV (log) RMSE: 0.141 ± 0.053

Tabla 17: Métricas finales

Métrica	Train log	Test log	Train orig	Test orig
RMSE	0.148	0.391	9.707	46.109
R²	0.739	0.202	0.797	0.013

Tabla 18: Bootstrap (coeficientes, 95 % CI)

Banda	Media β	IC 95% inf.	IC 95% sup.
log(B02/B07)	-0.5622	-0.6158	-0.5001

Ecuación (log)

log(sol_sus) = 4.2877 + -0.5625·log(B02/B07)

Ecuación (original)

sol_sus = 72.8024 · log(B02/B07)^-0.5625

Impacto del aumento del umbral ΔAIC en la selección de variables

Al aumentar ΔAIC a 4, el modelo se simplifica seleccionando solo la razón log(B02/B07), manteniendo un buen ajuste en entrenamiento (R² = 0.739). Sin embargo, el desempeño en la campaña de prueba es bajo, con R² cercano a cero, indicando poca capacidad de generalización. El análisis bootstrap muestra estabilidad en el coeficiente seleccionado.

Por lo tanto, un mayor ΔAIC favorece modelos más simples, pero no mejora la predicción en nuevas campañas, por lo que es necesario explorar modelos más complejos o incluir más variables.

2.9 Selección de bandas, razones y potencias para modelar sólidos suspendidos

En esta sección amplío el conjunto de variables predictoras al incluir no solo logaritmos simples de bandas espectrales, sino también sus potencias cuadráticas y cúbicas, así como razones logarítmicas entre bandas. El objetivo es capturar relaciones no lineales y combinadas que podrían mejorar la capacidad predictiva del modelo.

Utilizo una estrategia de selección hacia adelante basada en el criterio AIC (ΔAIC ≥ 4), priorizando modelos más parsimoniosos pero informativos. La validación se realiza con partición temporal (hold-out) y GroupKFold, y se evalúa la estabilidad del modelo mediante bootstrap de coeficientes (IC del 95 %). Finalmente, se comparan las predicciones con observaciones en escalas logarítmica y original.

Código

import numpy as np
import pandas as pd
from pathlib import Path
from sklearn.linear_model import LinearRegression
from sklearn.model_selection import GroupKFold
from sklearn.metrics import mean_squared_error, r2_score
import matplotlib.pyplot as plt
from IPython.display import display, Markdown

# -------------------------------------------------------------
# 1 · CONFIGURACIÓN
DATOS     = Path("datos/datos_sol_sus_acolite.csv")
DELTA_AIC = 4                # umbral de mejora
N_BOOT    = 2000
RNG_SEED  = 42

# -------------------------------------------------------------
# 2 · FUNCIONES
def aic_gauss(rss, n, k):              # k = coef + intercepto
    return n*np.log(rss/n) + 2*k

def fit_metrics(X, y):
    m = LinearRegression().fit(X, y)
    pred = m.predict(X)
    rss  = np.sum((y-pred)**2)
    n,p  = X.shape
    return {
        "AIC" : aic_gauss(rss, n, p+1),
        "RMSE": np.sqrt(mean_squared_error(y, pred)),
        "R²"  : r2_score(y, pred),
        "R²aj": 1 - (1-r2_score(y,pred))*(n-1)/(n-p-1),
        "mod" : m
    }

rmse = lambda a,b: np.sqrt(mean_squared_error(a,b))

# -------------------------------------------------------------
# 3 · CARGA Y PREPARACIÓN
df = pd.read_csv(DATOS)
df["fecha"] = pd.to_datetime(df["fecha"])

bandas = [c for c in df.columns if c.startswith("B")]
df = df.loc[(df["sol_sus"]>0) & (df[bandas]>0).all(1)].copy()

df["log_sol_sus"] = np.log(df["sol_sus"])

# 3.1  log(Bx)  +  (log(Bx))²,(log(Bx))³
log_cols   = []
poly_cols  = []
for b in bandas:
    log_b = f"log_{b}"
    df[log_b] = np.log(df[b])
    df[f"{log_b}^2"] = df[log_b]**2
    df[f"{log_b}^3"] = df[log_b]**3
    log_cols.append(log_b)
    poly_cols.extend([f"{log_b}^2", f"{log_b}^3"])

# 3.2  razones en log: log(Bi/Bj) = logBi - logBj
ratio_cols = []
for i,b1 in enumerate(bandas):
    for b2 in bandas[i+1:]:
        col = f"log({b1}/{b2})"
        df[col] = df[f"log_{b1}"] - df[f"log_{b2}"]
        ratio_cols.append(col)

cand_cols = log_cols + poly_cols + ratio_cols
print(f"Total de variables candidatas: {len(cand_cols)}")

X = df[cand_cols]
y = df["log_sol_sus"]

# -------------------------------------------------------------
# 4 · HOLD‑OUT TEMPORAL
last = df["fecha"].max()
m_te  = df["fecha"]==last
m_tr  = ~m_te
X_tr, y_tr = X[m_tr], y[m_tr]
X_te, y_te = X[m_te], y[m_te]
groups_tr  = df.loc[m_tr,"fecha"]

print(f"\nTRAIN campañas : {groups_tr.nunique()}")
print(f"TEST  campaña  : {last.date()}  ({m_te.sum()} muestras)")

# -------------------------------------------------------------
# 5 · FORWARD‑AIC (Δ ≥ {DELTA_AIC})
remaining = cand_cols.copy()
selected  = []

rss0   = np.sum((y_tr - y_tr.mean())**2)
best_aic = aic_gauss(rss0, len(y_tr), 1)

hist=[]; step=0
while remaining:
    scores=[]
    for col in remaining:
        cols = selected + [col]
        mtr  = fit_metrics(X_tr[cols].values, y_tr.values)
        scores.append((mtr["AIC"], col, mtr))
    scores.sort()
    new_aic, best_col, m_info = scores[0]

    if best_aic - new_aic >= DELTA_AIC:
        selected.append(best_col)
        remaining.remove(best_col)
        best_aic = new_aic
        step += 1
        hist.append({
            "Paso":step, "Variables":", ".join(selected),
            "AIC":best_aic, "RMSE":m_info["RMSE"],
            "R²":m_info["R²"], "R²aj":m_info["R²aj"]
        })
    else:
        break

if not selected:
    best_col = min(scores, key=lambda t:t[0])[1]
    selected = [best_col]
    hist.append({"Paso":1,"Variables":best_col,
                 **{k:v for k,v in scores[0][2].items() if k!="mod"}})

hist_df = pd.DataFrame(hist)

# Mostrar tabla historial con caption y sin índice
display(
    hist_df.style
    .format({"AIC": "{:.3f}", "RMSE": "{:.3f}", "R²": "{:.3f}", "R²aj": "{:.3f}"})
    .set_caption("Historial de selección")
    .hide(axis="index")
)

# Mostrar variables seleccionadas en negrita Markdown
display(Markdown(f"**Variables finales seleccionadas:** `{', '.join(selected)}`"))

# -------------------------------------------------------------
# 6 · VALIDACIÓN GROUP‑KFOLD
gkf = GroupKFold(n_splits=min(5,groups_tr.nunique()))
rm = []
for tr_idx, te_idx in gkf.split(X_tr[selected], y_tr, groups_tr):
    m = LinearRegression().fit(X_tr.iloc[tr_idx][selected], y_tr.iloc[tr_idx])
    rm.append(rmse(y_tr.iloc[te_idx], m.predict(X_tr.iloc[te_idx][selected])))

display(Markdown(f"**GroupKFold CV (log) RMSE:** {np.mean(rm):.3f} ± {np.std(rm):.3f}"))

# -------------------------------------------------------------
# 7 · MODELO FINAL Y MÉTRICAS
mod = LinearRegression().fit(X_tr[selected], y_tr)

pred_tr = mod.predict(X_tr[selected])
pred_te = mod.predict(X_te[selected])

metrics_df = pd.DataFrame({
    "Métrica"    : ["RMSE", "R²"],
    "Train log"  : [rmse(y_tr,pred_tr), r2_score(y_tr,pred_tr)],
    "Test log"   : [rmse(y_te,pred_te), r2_score(y_te,pred_te)],
    "Train orig" : [rmse(np.exp(y_tr), np.exp(pred_tr)), r2_score(np.exp(y_tr), np.exp(pred_tr))],
    "Test orig"  : [rmse(np.exp(y_te), np.exp(pred_te)), r2_score(np.exp(y_te), np.exp(pred_te))]
})

# Mostrar tabla métricas con caption y sin índice
display(
    metrics_df.style
    .format({
        "Train log": "{:.3f}", "Test log": "{:.3f}",
        "Train orig": "{:.3f}", "Test orig": "{:.3f}"
    })
    .set_caption("Métricas finales")
    .hide(axis="index")
)

# -------------------------------------------------------------
# 8 · BOOTSTRAP AGRUPADO (IC 95 %)
rng = np.random.default_rng(RNG_SEED)
dates = groups_tr.unique()
coef_bs=[]
for _ in range(N_BOOT):
    samp = rng.choice(dates, len(dates), replace=True)
    idx  = groups_tr.isin(samp)
    m_bs = LinearRegression().fit(X_tr.loc[idx,selected], y_tr.loc[idx])
    coef_bs.append(m_bs.coef_)
coef_bs = np.asarray(coef_bs)
ci_lo, ci_hi = np.percentile(coef_bs,[2.5,97.5],axis=0)

bs_df = pd.DataFrame({
    "Variable"    : selected,
    "Media β"     : coef_bs.mean(0),
    "IC 95% inf." : ci_lo,
    "IC 95% sup." : ci_hi
})

# Mostrar tabla bootstrap con caption y sin índice
display(
    bs_df.style
    .format({
        "Media β": "{:.4f}",
        "IC 95% inf.": "{:.4f}",
        "IC 95% sup.": "{:.4f}"
    })
    .set_caption("Bootstrap (coeficientes, 95 % CI)")
    .hide(axis="index")
)

# -------------------------------------------------------------
# 9 · ECUACIÓN DEL MODELO
alpha = np.exp(mod.intercept_)
expr_log  = " + ".join([f"{c:.4f}·{v}" for c,v in zip(mod.coef_,selected)])
expr_orig = " · ".join([f"{v}^{c:.4f}" for c,v in zip(mod.coef_,selected)])

display(Markdown("Ecuación (log)"))
display(Markdown(f"`log(sol_sus) = {mod.intercept_:.4f} + {expr_log}`"))

display(Markdown("Ecuación (original)"))
display(Markdown(f"`sol_sus = {alpha:.4f} · {expr_orig}`"))

# -------------------------------------------------------------
# 10 · PARIDAD TRAIN + TEST (log y mg/L) con estilo de colores y leyenda

plt.rcParams.update({"axes.titlesize":10,"axes.labelsize":8,
                     "xtick.labelsize":7,"ytick.labelsize":7})

TRAIN_CLR = "#9D50A6"
TEST_CLR = "#FFA24B"
ID_CLR = "#17A77E"

fig, ax = plt.subplots(1, 2, figsize=(9, 4))
# Escala log
ax[0].scatter(y_tr, pred_tr, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[0].scatter(y_te, pred_te, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims = [min(y_tr.min(), y_te.min()), max(y_tr.max(), y_te.max())]
ax[0].plot(lims, lims, "--", color=ID_CLR, lw=2.5, label="y = x")
ax[0].set(xlabel="log(sol_sus) obs.", ylabel="log(sol_sus) pred.", title="Escala log")
ax[0].legend(fontsize=7, frameon=False)
ax[0].grid(ls=":")

# Escala original
y_tr_lin, y_te_lin = np.exp(y_tr), np.exp(y_te)
pred_tr_lin, pred_te_lin = np.exp(pred_tr), np.exp(pred_te)

ax[1].scatter(y_tr_lin, pred_tr_lin, alpha=0.6, color=TRAIN_CLR, edgecolor="k", label="Train (●)")
ax[1].scatter(y_te_lin, pred_te_lin, alpha=0.6, color=TEST_CLR, marker="s", edgecolor="k", label="Test  (■)")
lims_o = [min(y_tr_lin.min(), y_te_lin.min(), pred_tr_lin.min(), pred_te_lin.min()),
          max(y_tr_lin.max(), y_te_lin.max(), pred_tr_lin.max(), pred_te_lin.max())]
ax[1].plot(lims_o, lims_o, "--", color=ID_CLR, lw=2.5)
ax[1].set(xlabel="sol_sus obs. (mg/L)", ylabel="sol_sus pred. (mg/L)", title="Escala original")
ax[1].legend(fontsize=7, frameon=False)
ax[1].grid(ls=":")

fig.suptitle(f"Predicho vs Observado · ΔAIC ≥ {DELTA_AIC}", fontsize=11, y=1.05)
fig.tight_layout(pad=1.2)
plt.show()

# Restaurar estilos para otras gráficas
plt.rcParams.update({"axes.titlesize":11,"axes.labelsize":9,
                     "xtick.labelsize":8,"ytick.labelsize":8})

Total de variables candidatas: 88

TRAIN campañas : 5
TEST  campaña  : 2025-06-09  (12 muestras)

Tabla 19: Historial de selección

Paso	Variables	AIC	RMSE	R²	R²aj
1	log(B02/B07)	-106.889	0.148	0.739	0.729
2	log(B02/B07), log_B8A^3	-114.331	0.126	0.812	0.797

Variables finales seleccionadas: log(B02/B07), log_B8A^3

GroupKFold CV (log) RMSE: 0.123 ± 0.034

Tabla 20: Métricas finales

Métrica	Train log	Test log	Train orig	Test orig
RMSE	0.126	0.244	8.384	33.133
R²	0.812	0.688	0.849	0.490

Tabla 21: Bootstrap (coeficientes, 95 % CI)

Variable	Media β	IC 95% inf.	IC 95% sup.
log(B02/B07)	-0.7848	-0.9849	-0.0678
log_B8A^3	-0.0036	-0.0075	0.0085

Ecuación (log)

log(sol_sus) = 4.1224 + -0.8370·log(B02/B07) + -0.0044·log_B8A^3

Ecuación (original)

sol_sus = 61.7089 · log(B02/B07)^-0.8370 · log_B8A^3^-0.0044

2.10 Selección de variables mediante Forward Selection con AIC usando validación cruzada

Tabla 22: Selección variables (AIC, RMSE, R², R²aj)

Variables	RMSE	R²	R²aj	AIC
['log_ratio_B04_B05']	0.184	0.782	0.747	-24.030
['log_ratio_B04_B05', 'log_B04']	0.161	0.829	0.764	-24.393

Variables finales: log(B04/B05), log(B04)

Ecuación en escala logarítmica:

\(\displaystyle \log(sol\_sus) = 4.459− 0.424\,\log\left(\frac{B04}{B05}\right) − 0.116\,\log(B04)\)

Ecuación en escala original:

\(\displaystyle sol\_sus = 86.402 \times \left(\frac{B04}{B05}\right)^{-0.424} \times B04^{-0.116}\)

Tabla 23: Desempeño Train y Test

Métrica	Train (log)	Test (log)	Train (mg/L)	Test (mg/L)
RMSE	0.123	0.231	12.177	12.847
R²	0.904	0.637	0.904	0.864
R² adj	0.897	0.547	0.897	0.829

Tabla 24: Intervalos de Confianza 95 % (Bootstrap)

Métrica	Train prom	Train 2.5%	Train 97.5%	Test prom	Test 2.5%	Test 97.5%
RMSE	0.115	0.083	0.144	0.171	0.062	0.244
R²	0.908	0.847	0.954	0.677	0.170	0.960
R² adj	0.901	0.836	0.951	0.596	-0.038	0.950

2.11 Selección de variables mediante Forward Selection con AIC usando bootstrap

Tabla 25: Selección variables (Bootstrap completo)

Variables	RMSE	R²	R²aj	AIC
['log_ratio_B04_B05']	0.174	0.813	0.808	-136.798
['log_ratio_B04_B05', 'log_B02']	0.146	0.867	0.860	-148.468

Variables finales: B04/B05, B02

Ecuación en escala logarítmica:

\(\displaystyle \log(sol\_sus) = 4.459−0.371\,\log(B04/B05)−0.091\,\log(B02)\)

Ecuación en escala original:

\(\displaystyle sol\_sus = 86.402\left(\frac{B04}{B05}\right)^{-0.371} \times B02^{-0.091}\)

Tabla 26: Desempeño Train y Test

Métrica	Train (log)	Test (log)	Train (mg/L)	Test (mg/L)
RMSE	0.129	0.206	12.528	11.233
R²	0.894	0.712	0.899	0.896
R² adj	0.886	0.640	0.891	0.870

Tabla 27: Intervalos de Confianza 95 % (Bootstrap)

Métrica	Train prom	Train 2.5%	Train 97.5%	Test prom	Test 2.5%	Test 97.5%
RMSE	0.120	0.089	0.147	0.161	0.065	0.221
R²	0.897	0.821	0.945	0.713	0.191	0.950
R² adj	0.889	0.808	0.941	0.641	-0.011	0.937

2.12 Procesamiento de datos geográficos

2.12.1 Modelo gkf

\(\displaystyle \text{sol\_sus} = 61.7089 \cdot \left(\frac{B02}{B07}\right)^{-0.8370} \cdot \left(B8A^3\right)^{-0.0044}\)

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 76.45%

--- Umbral: 0 ---
Estadísticas del índice de agua:
  Min: -12.4261, Max: 17.4322, Media: 0.6814
Estadísticas de la máscara de agua:
  Valores únicos: [0 1]  (1=agua, 0=tierra, 255=nodata)
Estadísticas del parámetro sólidos suspendidos:
  Min: 0.1966, Max: 12618.6895, Media: 74.7291

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 76.45%
Umbral sugerido: 0.330

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.87%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.74%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.73%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.71%

2.12.2 Modelo vc

\(\displaystyle \mathrm{sol\_sus} = 86.402 \cdot \left(\tfrac{B04}{B05}\right)^{-0.424} \cdot B04^{-0.116}\)

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 76.45%

--- Umbral: 0 ---
Estadísticas del índice de agua:
  Min: -12.4261, Max: 17.4322, Media: 0.6814
Estadísticas de la máscara de agua:
  Valores únicos: [0 1]  (1=agua, 0=tierra, 255=nodata)
Estadísticas del parámetro sólidos suspendidos:
  Min: 16.0556, Max: 2994.8086, Media: 113.4642

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.87%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.74%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.73%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.71%

2.12.3 Modelo bt

\(\displaystyle \mathrm{sol\_sus} = 86.402 \cdot \left(\frac{B04}{B05}\right)^{-0.371} \cdot B02^{-0.091}\)

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 76.45%

--- Umbral: 0 ---
Estadísticas del índice de agua:
  Min: -12.4261, Max: 17.4322, Media: 0.6814
Estadísticas de la máscara de agua:
  Valores únicos: [0 1]  (1=agua, 0=tierra, 255=nodata)
Estadísticas del parámetro sólidos suspendidos:
  Min: 19.9269, Max: 51112.1016, Media: 139.5654

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.87%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.74%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.73%

Índice de agua → min: -12.426, max: 17.432
% de píxeles clasificados como agua: 74.71%

Referencias

[1]

D. C. R. Ramírez, «Método de Estimación de Sólidos Suspendidos Totales como Indicador de la Calidad del Agua Mediante Imágenes Satelitales». Universidad Nacional de Colombia, Facultad de Ciencias Agrarias, 2017.

[2]

G. D. J. L., J. Stephan, y Delance-Martinic., «Determinación del parámetro sólidos suspendidos totales (sst) mediante imágenes de sensores ópticos en un tramo de la cuenca media del río Bogotá (Colombia).», UD y la Geomática, pp. 19-27, 2014, doi: 10.14483/udistrital.jour.udgeo.2014.9.a02.

[3]

A. Cruz-Retana, C. Fonseca-Ortiz, R. Becerril-Piña, M. A. Gómez-Albores, M. Hernández-Téllez, y R. Arévalo-Mejia, «Characterization of spectral reflectance and TSS concentration in continental water bodies worldwide», vol. 1. pp. 4-18, 2023.

[4]

E. E. C. Vargas, «Modelamiento de Relaciones entre Parámetros Fisicoquímicos y Microbiológicos en Aguas de la Bahía Interior del Lago Titicaca-Puno (Perú) mediante Árboles de Predicción», Revista Tecnica De La Facultad De Ingenieria Universidad Del Zulia, vol. 44, pp. 154-168, ago. 2021, doi: 10.22209/rt.v44n3a02.

[5]

M. Moeini, A. Shojaeizadeh, y M. Geza, «Supervised machine learning for estimation of total suspended solids in urban watersheds», Water (Switzerland), vol. 13, ene. 2021, doi: 10.3390/w13020147.

[6]

L. S. Kupssinskü, T. T. Guimarães, E. M. D. Souza, y D. C. Zanotta, «A method for chlorophyll-a and suspended solids prediction through remote sensing and machine learning», Sensors (Switzerland), vol. 20, abr. 2020, doi: 10.3390/s20072125.

Notas

Aguas lénticas.↩︎
d = prueba estadística de Durbin-Watson.↩︎